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SUMMARY 


Development of a flow prediction method for rocket turbopumps is discussed in this report. A 
description is given of the complex nature of the flowfield existing in turbopumps. Examples 
are given to illustrate that both physics based models and analytical calculation procedures based on 
Computational Fluid Dynamics (CFD) have resulted in progressive advancements of design 
procedures used in turbopumps. Limitations of the state— of— the— art design procedures are outlined 
for which the present work was conducted to significantly enhance the design methodology. 

A CFD code developed at NASA ARC was used as the base code in the present work. In its initial 
form the CFD code could compute unsteady turbulent flow through an axial flow turbine stage. 
Governing equations and numerical procedure used in the CFD code is documented in detail. The 
turbulence model in the code was modified to facilitate computation of transitional flows and to 
ac co unt for extra rates of strain, such as rotation, three — dimensionality, surface curvature and surface 
roughness. Boundary conditions in the code were modified to facilitate computation of surface heat 
transfer coefficients and to allow computation through multistage turbomachines. The code was 
modified to permit simulations of flow through airfoil rows with flowpath convergence and divergence 
and to include radial flow turbomachines. 

Extensive work was conducted to demonstrate that the CFD code yields good estimates of airfoil 
loadings, heat transfer coefficients, boundary layer parameters, losses, endwall secondary flows and 
tip leakage flows. Benchmark quality data, obtained from two— and three-dimensional cascades, are 
used in the code verification process. The ability of the base code to compute time-averaged and 
unsteady flow through a turbine stage had been previously demonstrated by the originator of the code 
at NASA ARC. Additional computations were conducted by NASA ARC personnel concurrent to the 
present program to demonstrate that the unsteady and time-averaged flow prediction capabilities 
of the code could be improved by utilizing more refined grids and by accounting for a more realistic 
airf oil count in axial flow turbines. Tb avoid duplication of work and to build on the enhanced strength 
of the code, work was conducted to demonstrate that the present code yields a more realistic estimate 
of unsteady loads in turbines than unsteady Euler codes; the latter codes are currently being used in 
the design procedures for rocket turbopumps. 

Computations were conducted to demonstrate that the modified code with improved turbulence 
models, developed in the present program, yields realistic estimates of unsteady and time -averaged 
losses in multi-stage turbomachines. Work conducted in the present program and other simulations 
conducted by the authors during the program period indicated that the present code, operated in a 
two-dimensional mode (but modified to account for stream-tube variation effects), is a cost 
effective alternative to full three-dimensional calculations. This approach permits realistic 
predictions of unsteady loadings and losses for multistage machines, allowing design engineers more 
time for design optimization studies. The predictive capabilities of the present code were 
demonstrated by computing flow through a radial impeller and a multistage axial flow turbine. 

The work conducted in the present program was supported by NASA MSFC Contract 
#NAS8— 36950. 
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1. INTRODUCTION 


TUrbopumps for future rocket engines will require operating lives and performance levels well 
above those of present day units. Due to the very hostile environment that exists during the actual 
operation of a turbopump, reliable data that would permit identification of specific problems and thus 
guide future designs, is extremely difficult to obtain. In recognition of this situation, NASA has 
sponsored Computational Fluid Dynamics (CFD) computer code development programs which, 
when verified against basic data, may be used to identify ways to improve the current turbopumps and 
provide a basic understanding of the flow phenomena which could lead to better designs of future 
turbopumps for rocket engine applications. This report presents results from one of these programs, 
funded by NASA MSFC under contract # NAS8- 36950 on the development, modification, 
verification and application of a CFD code for turbopump design application. 
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2. BACKGROUND 


Although real flows in rocket turbopumps are highly unsteady, steady flow calculation methods 
are used during the design process. The unsteadiness is induced by the temporal and spatial variation 
of inlet temperatures and by aerodynamic interaction between rotating and stationary components. 
In the design process the effects of unsteadiness, at best, are accounted for through empirical 
correlations derived from available experimental data. Lack of physics in these empirical correlations 
invariably yield non— optimal designs requiring extensive efforts which add time and cost to the 
development process. Examples are given below to demonstrate that in addition to the fluctuations 
in the aerodynamic loads, unsteadiness affects time averaged losses and heat loads which must be 
accounted for in the design process. 

The interaction which occurs in the SSME turbopump between the rotating pump impeller and 
the stator diffuser vane and also in the pump turbine between the rotor and the stator can generate 
unsteady forces through various mechanisms which occur at widely different time and length scales. 
Some of these mechanisms are summarized below (Greitzer (1987)) 



Unsteadiness Due to: 

i) Potential field interaction due to 
relative motion of airfoils 

ii) Wakes 

iii) Large scale vortical organized 
flow structures 

iv) Tfcmperature distortions 

v) 

vi) 

vii) 


chord 

5 X lO" 5 sec 

chord 

5 X lO' 5 sec 

chord/span 

5 X 10” 5 sec 

gap/span 

5 X 10” 5 sec 


Tbrbulence chord 

Rotating stall circumference 5 X 10“ 1 sec 

Surge length 5 X 10 * sec 

Typical row assumed to have 60 airfoils at gap/chord of one and operating at 36000 RPM 


These mechanisms are governed by different types of fluid dynamic phenomena. In order for any 
CFD code to capture these unsteady interactions, the responsible fluid dynamic phenomena tabulated 
above must be accurately modeled. A brief discussion on each of these phenomena, except rotating 
ctaii and surge, is given below. Rotating stall and surge typically occur in high pressure ratio 
compressors, and they may not contribute to the unsteadiness prevalent in the SSME turbopump. 
Unsteadiness due to turbulence is primarily embedded in the unsteadiness generated by wakes and 
large scale vortical organized structures; however, at the first row of airfoils, it can primarily be treated 
as a steady— state phenomenon using turbulence models in Reynolds— Averaged Navier— Stokes 
codes. 
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2.1 Potential Flowfield Interaction 


The potential flowfield interaction in the turbopump components arises because the vane or the 
impeller is subject to a time— varying pressure field influence resulting from relative movement of 
adjacent airfoil rows. This time— varying force has a time scale on the order of the blade passing 
frequency and could cause structural fatigue in and of itself. In addition, the time— varying pressure 
field is capable of creating a transient flow disturbance, such as separation, in the vane or impeller 
passage; a possible interactive feedback mechanism could sustain or amplify the unfavorable 
interaction between the rotating and stationary components. 


The potential interaction induced by relative motion of the adjacent airfoil rows is shown in 
Figure 2.1.1. This figure indicates that the pressure field (waves) due to the potential around an airfoil 
row extends both upstream and downstream of the airfoil. Typically the strength of the field decays 
over a length scale equal to the pitch/chord of the cascade. Unsteadiness effects due to relative motion 
of the pressure field, in both the upstream and the downstream rows increases with decreasing axial 
distance between adjacent rows. 



Figure 21.1 Potential Flow Pressure Gradients. 


The magnitude of loading variation for a turbine stator and rotor at mean radius as measured 
in the United Tfechnologies Research Center (UTRC) Large Scale Rotating Rig (LSRR) by Dring et 
al. (1982) is illustrated in Figure 2.1.2. The unsteady pressure variation on the stator caused by the 
upstream potential influence of the downstream rotor was as much as 15% of the steady stator exit 
dynamic pressure (Figure 2.1.2(a)). The pressure variations on the downstream rotor were measured 
when it was located at 15% axial chord downstream of the stator. The amplitudes of the pressure 
si gnals for the rotor were of the order of 80% of the relative steady state dynamic pressure. Data were 
also acquired when the gap between the upstream stator and downstream rotor was of the order of 
60% chord of the upstream stator. The amplitudes of the unsteady pressure signals on the rotor at 
hi gh er axial gap condition were measured to be about one— half those at the lower axial gap. The 
diff erence in the unsteady pressure amplitude on the rotor for the two axial gaps indicates a decay rate 
of pressure amplitude with distance which is far less than the rate of decay of potential influence. This 
implies that a large portion of the unsteadiness in pressure on the downstream rotor airfoil is due to 
wakes from the upstream rotor. A discussion on the physics of wakes as it affects downstream airfoil 
rows is given in the following subsection. 



Figure 21.2 Range of Instantaneous Pressure Distributions Measured in the 

United Technologies Research Center Large Scale Rotating Rig 
(UTRC LSRR) on Stator and Rotor Airfoils. 

The impact of potential flow interaction in a turbopump was numerically demonstrated by 
Rangwala et al. (1991) for the Generic Gas Generator (G3) turbine. Calculations were conducted for 
the mean section of the G3 turbine first stage, for three axial gaps between the first stator and the rotor, 
by using a two dimensional (2D) unsteady Navier— Stokes code. Results from these calculations 
(Fig.2.1 3) clearly show that time-averaged diffusion on the upstream stator is a function of the axial 
gap between the vane and the blade. The upstream vane is, therefore, affected by the downstream 
blade through the potential effect. The above calculations also showed that the turbine efficiency 
could be improved by up to 1% by operating the turbine at the larger gap. Although experimental data 
verifying these numerical results are not yet available, indications are that there is an optimum gap 
between adjacent airfoil rows. 



0.0 Normalized Axial Distance 1.0 0.0 Normalized Axial Dlatance 1.0 


Figure 21.3 Simulations (Rangwala et aL (1991)) Conducted for a 
TUrbine Stage by Using a 2D Unsteady Navier- Stokes 
Code Show That Time-Averaged Diffusion on the 
Upstream Vane is Influenced by the Axial Gap Between 
the Vane and the Rotor. 
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22 Wake & Temperature Streak Interaction 

The circumferential variations in the velocity field downstream of the first stators for turbines 
are normally generated by the drag on the airfoil and endwall surfaces which causes a reduction in 
velocity and increases in the turbulence levels in the low velocity regions. In some flow situations, 
especially for airfoil rows downstream of a combustor, high velocity jets exist due to large 
circumferential gradients in temperatures. The effects of these upstream velocity variations can be 
simply illustrated through the use of velocity triangles (Butler et al. (1986)), as shown in Figure 2.2.1 
for a turbine. This figure shows that the lower velocity fluid has a normal velocity component towards 
the suction side of the downstream airfoil indicating that the high turbulence, low momentum fluid 
from the upstream airfoil wake will migrate towards the suction side of the airfoil. In a similar manner, 
high velocity (high temperature) fluid will migrate towards the pressure side of the downstream 
airfoil. This preferential migration of fluid particles has three effects: 

i) Alterations in the boundary layer characteristics of the airfoil through its effect on the 
transition process. 

ii) Variation in the secondary flow generation for downstream passages. 

iii) Redistribution of stagnation enthalpy. 



U MMOL WOO 


Figure 221 Rotor Inlet Gas Temperature Distortion Causes Large 

Variation in Rotor Airfoil Incidence Angle. Simple 
Calculations Conducted for Hot to Cold Temperature 
Ratio of 1.7 Indicates Incidence Angle Variation of 12* 
and 40* for Typical High and Low Flow Coefficient. 
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2*2.1 Effect of Upstream Wakes on Losses and Heat Loads 

Detailed experimental investigations (Hodson (1983), Pfeil & Herbst (1979) Doorley & 
Oldfield (1984), Dring et a 1. (1982)) have been conducted to identify the influence of upstream 
periodic wakes on boundary layer characteristics. The characteristic of boundary layers on airfoil 
suction sides, affected by periodic movement of the upstream wakes, was clearly illustrated by Doorley 
& Oldfield (1984). Time resolved heat transfer data were acquired in this investigation in a stationary 
cascade at a low and a high background turbulence level with and without an upstream rotating rod 
that simulated wakes from upstream airfoils. These data, plotted in Figure 2.2.2(a), indicate that 
upstream wakes have relatively little effect on turbulent boundary layers and significant effects on 
laminar boundary layers. This figure shows that time— averaged heat transfer and losses for turbine 
airfoils, which normally have large regions of laminar flow in a steady flow environment, are expected 
to increase in an unsteady environment. 


•I WITHOUT UNSTEADINESS W WITH UNSTEADINESS 



\ I 

i*~VANE PASSING-^ 
I TIME I 

i 


TIME 

<Q— *- DATA WITH LOW INLET TURBULENCE LEVELS 
<g)— »- DATA WITH HIGH INLET TURBULENCE LEVELS 

Figure 112(a) Measured Time-Resolved Heat Transfer on a TUrbine Airfoil 
Suction Side (Doorfy et aL (1984)) at Two Background 
TUrbulence Levels in an Unsteady Environment Shows a Larger 
Effect on a Laminar Boundary Layer and Little Effect on a 
TUrbulent Boundary Layer. 

The time— averaged effect of upstream wakes on the boundary layer thickness and heat transfer 
coefficient for the suction sides of two separate rotor airfoils (Hodson (1983), Blair et al. (1988), 
Sharma et al. (1988)) are shown in Figure 2.2.2(b). Also shown in this figure are the data obtained for 
those airfoils in steady cascade configurations and calculated values from a boundary layer code. The 
steady cascade data are shown to yield good agreement with transitional calculations. The 
time— averaged data, however, lies between the transitional and fully— turbulent calculations. This 
figure indicates that the nature of transition is influenced by the periodic variation of turbulence 
imposed by wakes from upstream airfoil rows. In addition, losses and heat loads in an unsteady 
environment are larger than those measured in steady cascade configurations. For typical turbine 
airfoils, losses induced by unsteady effects may be on the order of 25% to 100% of the losses for those 
airfoils in a steady flow environment (Sharma et al. (1988)). 
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Figure 112(b) Measured Streamwise Distribution of Time-Averaged Stanton 
Number of Blair et aL (1988) and Sharma et aL (1988) and 
Boundary Layer Thickness ofHodson (1983) Show Larger Values 
in an Unsteady Environment Than in a Steady Cascade 

Configuration. 

233 Effect of Upstream Wakes on Secondary Flows 

In addition to affecting the characteristics of airfoil boundary layers, wakes from upstream airfoil 
rows also affect the generation of secondary flows (Sharma et al. (1983,88,90) as discussed below. 

Unsteady experimental data for a rotor passage, obtained by using high response probes in the 
UTRC LSRR (Sharma et al. (1983)) are shown in Figure 2.23. This figure shows three instantaneous 
contour plots of relative total pressure coefficient upstream and downstream of the rotor passage. 
Large variations in the exit flow structures are seen in the figure for the three different inlet conditions. 
The exit flow field (Figure 2.2.3(a)) shows three distinct vortices due to the hub and tip secondary flows 
and the tip leakage effects. Without the tip leakage vortex, the flow field shown in Figure 2.23(a) is 
similar to the one expected for this airfoil in a steady cascade environment. The tip leakage vortex for 
the rotor shows least variation (Figures 233(a), (b), and (c)) indicating that the leakage phenomenon 
is not influenced by upstream circumferential distortions. The hub secondary flow vortex shows the 
largest variation transforming from a distinct structure in Figure 233(a) to a diffused structure in 
Figure 23.3(b), and becoming almost non-existent in Figure 2.23(c). This indicates that the 
secondary flow generation mechanisms, especially at the hub, are strongly influenced by the upstream 
circumferential distortions such as wakes. The overall variation in the size of the tip secondary flow 
vortex is smaller than that of the hub vortex but larger than the leakage vortex. 
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Figure 113 Secondary Flow Structures Downstream of a Rotor (Sharma et aL (1988)) 

Obtained From Unsteady Measurements Show Large Variation in Their 
Size, Indicating Effects of Upstream Stator Wakes. 

The periodic variation in the size and strength of the secondary flow vortices observed in this 
experimental investigation shows almost 40% variation in the secondary flow losses for the rotor 
passage. 

2£3 Eff ects of Upstream Temperature Streaks on Segregation of Hot and Cold Air in 
Tkirbine Rotors 

Results from an experimental investigation, conducted to quantify the influence of burner 
hot streaks on segregation of hot and cold air in turbine rotors, are discussed below. 

In this investigation, experimental data were acquired in the UTRC LSRR by introducing 
temperature streaks upstream of the first stator (Figure 2.2.4(a)). TWo types of temperature profiles 
were generated upstream of the first stator, these being: 
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Figure 2.2.4(a) Schematics of the Experimental Apparatus Used 
to Simulate the Redistribution of Hot Streak in 
the TUrbine Rotor (Butler et oL (1986)). 


i) Hot streak with a circular cross-section to yield a temperature profile both in the radial and 
the circumferential direction, some of the results from this investigation were reported by 
Butler et al. (1986). 

ii) Hot streak with a rectangular cross-section to yield a radially uniform profile that had 
temperature gradients mainly in the circumferential direction (Sharma et al. (1990)). 

The hot air in these experiments was seeded with Carbon Dioxide (CO 2 ) to facilitate 
measurements of its migration in the turbine by using a gas sampling technique. The temperature 
patterns at the exit of the first stator for these tests are given in Figure 2.2.4(b) which indicates 
relatively small mixing in the first stator. This result is expected being compatible with the Munk and 
Prim (1947) principle. Spanwise distributions of axisymmetric CO 2 concentration profiles at inlet to 
the rotor,measured by using a rotating probe,for these two tests are given in Figure 2.2.4(c). The 
circular streak generate a parabolic concentration (temperature) profile, whereas the rectangular 
streak generated a radially uniform profile. 



Figure 114(b) Contour Plots of Normalized CO 2 Concentration 
Downstream of the First Stator in the UTRC LSRR 
Obtained with Circular and Rectangular Hot Streaks; 
High Values Imply High Temperatures. 


10 



NORMALIZED SPAN 

Figure 114(c) Spanwise Distribution of Normalized CO 2 Concentration 
Profiles (Indicators ofTemperatures) Measured in the Rotor 
Frame for the Circular and Rectangular Hot Streaks. 

Measured concentration of CO 2 on the rotor airfoil surfaces,for the rectangular hot streak, are 
shown in Figure 2.2.5. This figure indicates higher levels of CO 2 concentration on the rotor airfoil 
pressure side than on the suction side. Similar results were obtained with the circular hot streak as 
riis n ic<yrf by Butler et al. (1986). These results demonstrate that pressure sides of rotor airfoils 
operate at substantially higher temperatures than the suction sides when the incoming flow has a 
circumferentially non-uniform temperature profile. 



Figure 115 Larger Time-Averaged CO 2 Concentration (Temperatures) 

Measured on the Pressure Side of the Rotor Airfoil Relative to 
the Suction Side Indicate Segregation of Hot and Cold Air. 

Interpretation of data obtained from an aircraft engine environment indicates that the pressure 
sides of first rotor airfoils in a high pressure turbine can operate anywhere between 
100-700* F hotter than suction sides. These temperature differences between the two sides of airfoils 
can cause significant durability problems for airfoils. Large amounts of cooling air are required to 
accommodate these temperature levels resulting in reduced efficiency of the cycle and increased 
specific fuel consumption of aircraft engines. 
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23 Interaction Due to Large Scale Organized Vortical Flow Structures 

The real flow in the turbopump turbines consists of large scale organized flow structures formed 
by viscous flows in large deflection airfoil rows. Although these flow structures may affect the unsteady 
flow fields in turbines more than either the potential or wake interaction, this phenomenon has not 
been discussed in the open literature until recently by Sharma ct al. (1988,1990). Shortage of available 
data and complexity of the problem it introduces in the analytical treatment have contributed to this 
gap. Ihrbine cascade flow visualization experiments conducted in a water tunnel at University of 
Connecticut illustrate the organized structures formed in a simple configuration as shown in Figure 
23.1. The flow structures formed in the cascade were visualized by utilizing hydrogen bubbles and 
laser li ghting techniques. The deformation of two parallel horizontal lines of bubbles at the inlet and 
exit planes of the cascade is shown in the above figure. This figure clearly shows how the fluid particles 
contained between the two lines deform into vortical structures constituting almost 25% of the total 
airfoil passage. In a three dimensional unsteady flow field, these vortices should have a large influence 
on turbine performance, durability and structural integrity. 



Figure 2.3.1 Two Parallel Horizontal Lines Upstream of the Cascade Distort 
Into Vortical Structures at the Leading and Trailing Edges of the 
Cascade. Flow Visualization Tests Conducted at U. of Connecticut 
by Pratt & Whitney. 
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The formation and development of these vortices in turbine cascade passages were extensively 
investigated at UTRC (Langston et al. (1977)) in a large scale cascade rig where detailed aerodynamic 
performance, loadings and external heat loads were measured. The effect of these flows^rtweson 
the airfoil loading is shown in Figure 2.3.2; data were obtained on the cascade airfoil for three different 
inlet boundary layer profiles resulting in three different magnitudes of secondary flow vortices. Large 
spanwise gradients in airfoil loadings are generated with an increase in the incoming boundary layer 
thickness. This also results in an increase of the size of the vortex. The variation m the pitch averaged 
gas angle profile at the exit of the cascade increases from 5 to 25 degrees for the two extreme inlet 
profiles. Since these angle variations are introduced by vortical motions, they would cause similar 
variation in the circumferential direction, thus causing unsteady incidence variation on the following 
row of airfoils. The circumferential distortions introduced by these vortices are two to three times 
larger than those introduced by 2D wakes. 
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Figure 13.2 Change in Airfoil Loadings and Exit Gas Angles as Affected by Cascade 

Inlet Boundary Layers Langston’s Data (Langston et aL (1977), 
SharmaetaL (1990)). 
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Figure 2.3.2 Measured Streamline Patterns on the Airfoil Suction Side Indicating 

the Change in the Size of Secondary Flow Vortices Due to Inlet 
Boundary Layer. 

In addition to loading and exit angle variations the presence of secondary flow vortices also 
causes substantial effects on airfoil external heat loads as shown in Figure 233. This figure shows 
measured heat transfer coefficients on the midspan of the airfoil for thin and nominal inlet boundary 
layers; two different magnitudes of secondary flow vortices are represented. Although the airfoil 
loading at the midspan region for the two test conditions is the same, the heat load for the cascade 
with the larger secondary flow vortex is lower by 70%, indicating that the vortex causes substantial 
alterations to the structure of turbulence (Sharma & Graziani (1982)). The physical mechanisms 
governing the generation of secondary flow structures in steady cascade flows are well established 
(Langston et al. (1977) and Sharma & Butler (1986)), as depicted in Figure 23.4. The mechanisms 
governing the generation, development and transport of these vortices in unsteady turbomachinery 
environments, however, still require further work, especially from a predictive point of view. 
Extensive unsteady data acquired in the UTRC LSRR highlight important features of the secondary 
flow generation in the multistage turbine environment. These features have important implications 
on the development of analytical predictive models. 




Figure 23.3 Endwall Secondary Flow Vortex Affects External Heat Loads on 

Airfoil Suction Surface at Mid-Span (Sharma and Graziani 
(1982) ). Midspan Loadings in the 71 vo Tests are the Same. 



Figure 23.4 Cascade Endwall Flow Structure (Sharma &. Butler (1 986)) 

An e xam ple of large scale structures from the rotor influencing the flows in downstream stator 
airfoil rows is given in Figure 2.3.5. This figure shows the total pressure field downstream of a second 
stator which represents variations in the flowfield as one full rotor passage translates over one full 
second stator passage. The structures of the wakes and vortices of the second stator are only slightly 
affected by the upstream flow distortions, in contrast to information in Figure 2.2.3. Here, upstream 
flow distortions strongly influences the flow structure downstream of the rotor. A posable 
interpretation of the data shown in Figures 2.2.3 and 23.5 is that the effect of upstream flow 
distortions is enhanced by the rotation effects. The mid-passage flow region of the second stator, 
however, is strongly influenced by the upstream flow distortions. At certain time locations the vortices 
from the upstream rotor appear downstream of the second stator, completely unaffected by the stator. 
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Although these vortices contain fluid with large magnitudes of turbulence intensity and stresses, they 
show little evidence of mixing and thus indicate that these would have a large impact on the mid- span 
flow field of the second stator. 
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Figure 23.5 Unsteady Instantaneous Total Pressure Loss Coefficient Downstream of 
the Second Stator Indicate that Rotor Secondary Flow Vortice Periodically 
Persist Through the Second Stator (Sharma and Syed (1991)). 

Measured time averaged heat transfer data obtained at the midspan location of this airfoil is 
shown in Figure 23.6, together with theoretical predictions that were based on measured midspan 
airfoil loadings. The streamwise gradients of the measured and predicted heat transfer coefficients 
on the airfoil suction side differ from each other in the turbulent region; influence of the secondary 
flow vortices convecting through the midspan region of the stator on the airfoil boundary layer is 
apparent. 
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Figure 2.3 6 Analytical Predictions Underestimate External Heat Load on the 

Second Stator Airfoil Suction Surf ace for Blair et aL (1989) tests. 

This Result is Opposite to the Experiments and Predictions in 
Steady Cascade Configurations in Figure 2.2.2(b) and 2.3.3. 

The discussion presented in this section dearly demonstrates that large scale organized 
structures constitute a significant portion of the overall flow field in multistage turbomachines. The 
inddence angles and unsteadiness induced by these structures are larger than those induced by wakes 
and potential interactions for typical rocket turbines. No model or discussion of die existence of these 
flow structures and their influence on the unsteady flow field have appeared in literature, mdica g 
the complexity of the problem and infancy of the models. The development of a reliable P r fo ,ctl 
method, P suchas that discussed in the present program, will provide a valuable tool to turbopump 
designers, permitting optimization of the turbopump for improved durability, structural integrity and 

performance. 

In summary, unsteadiness effects in turbomachinery are caused not only by potential flow and 
two dimensional wake interactions but also by large scale organized structure and tem^rature 
distortions. Of these phenomena, large scale organized structures have been given little attention in 
the overall turbomachinery design and analysis procedures. These structures can cause higher 
unsteadiness in the airfoil rows than either potential flows or wakes. 

A review of the hierarchy of CFD codes used in the turbomachinery design and analyses modes 
is given in the following subsection, indicating how the work conducted under the present program 
permits improved analysis of turbopump flowfields to yield designs which are more efficient, durable 

and structurally sound. 


17 




14 Application of CFD Codes in Itirbomachinery 

Significant improvements in durability, structural integrity and performance of turbomachines 
have been realized over the past twenty years. One of the major contributors to these achievements 
has been the successful close coupling of analytical models describing the salient physical mechanism 
with state-of-the-art Computational Fluid Dynamics (CFD) codes. A brief review of the design 
process evolution with special emphasis on the use of CFD codes is given in the following subsection 
and schematically depicted in Figure 2.4.1. This will highlight how the work conducted under the 
present program has contributed to the enhancement of the current design process which will lead to 
improvements in durability, structural integrity and performance of rocket engine turbopump units. 



Figure 24.1 Application of Computational Fluid Dynamics Codes Have Resulted in 

Improved Performance. Further Performance Improvements Are 
Possible With Unsteady Code Applications. 

2.4.1 State -of- the- Art 

The actual flow fields in turbomachines are highly complex. They consist of laminar, transitional 
and turbulent boundary layers on airfoil surfaces and secondary flow vortices in the endwall regions. 
These vortices are formed by the incoming total pressure distortions in endwall boundary layers and 
tip leakage flows. These complex flow patterns are strongly influenced by three-dimensional 
pressure fields within the airfoil channels and relative movement of adjacent airfoils rows, as well as 
incoming time varying temperature patterns. The resultant flow field exhibits strongly viscous and 
highly time dependent characteristics, which necessarily create in transient thermal and aerodynamic 
loads on airfoil surfaces. 

In the early stages of turbomachinery design, before the advent of sophisticated computer 
systems, engineers relied on one -dimensional concepts and simple correlations derived from 


18 





extensive experimental data to account for loss generation mechanisms. These correlations, albeit 
relatively crude, accounted for the various flow field characteristics in a global sense and allowed 
engineers to generate consistent designs. A large volume of data acquired during the process of 
turbomachinery development was used to establish design criteria for airfoil loadings and radial 
distribution of flow, as shown in Figure 2.4.2. This one- dimensional approach to turbomachinery 
design often resulted in expensive and time-consuming product development. 




Figure 14.2 Low Loss Airfoil Design Criteria Established Through Extensive 
Data Base Review. No Leading Edge Diffusion, Large Acceleration 
Regions, Small Diffusion Regions. 


2.4.2 First Generation of CFD Code for Ttarbomachinery 


With the availability of computers in the mid 1960’s, CFD codes were developed to assist in the 
turbomachinery design process. The first generation of CFD codes used in the design solved 
two-dimensional equations both in the blade— to— blade (Caspar et al. (1979), Denton (1975) Ni 
(1982)) and the radial direction (Novak & Hearsay (1976)). Initially, these codes were used to analyze 
flowfield to ensure that the design intention was achieved with reduced hardware testing. TWo 
dimensional boundary layer codes were subsequently developed (Crawford & Kays (1976), 
McDonald & Fish (1972), Patankar & Spalding (1970) that permitted external heat load and 
aerodynamic loss calculations. Benchmark quality experiments (Blair & Werle (1980,81), Blair 
(1982), Sharma et al. (1982)), identifying the basic physics of boundary layers for typical turbine 
airfoils, were also conducted to develop improved turbulence models needed for the boundary layer 
codes. 


An example of results from these experiments is shown in Figure 2.4.3, where the effect of 
mains tream pressure gradient on transitional boundary layers is quantified. A simple algebraic 
turbulence model developed on the basis of these data yielded a much better estimate of the external 
heat loads on turbine airfoils than the more complex K— E turbulence models shown in Figure 2.4.4. 
The 2D inviscid and boundary layer codes can then be utilized to optimize shape of an airfoil during 
design, as shown in Figure 2.4.5. 
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Figure 24.3 Mainstream Pressure Gradient Effects More Pronounced on 
Thermal Boundary Layers Than on Momentum Boundary 
Layers in Transitional Region. 
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Figure 24.4 Transition Model of Sharma (1986) Gives Better Estimate of 
Heat Load for Consigny and Richard’s Airfoil (1982) Than the 
Launder— Jones (1972) Two Equation Model 



Figure 24.5 2-D Viscous and Inviscid Codes Used to Optimize 

Shapes of Airfoil for Improved Performance. 


The use of boundary layer calculation methods and basic experimental data also allow 
qualitative estimates of the effect of unsteadiness on turbine airfoil boundary layer development. 
Several specific experiments (Dring et al. (1982), Pfeil and his co-workers (1979,82), Hodson (1983), 
Sharma et al. (1983,88,90), Langston et al. (1977), Joslyn et al. (1982)) contributed to the development 
of the unsteady loss prediction models. Results from one of these experiments shown in Figure 2.4.6 
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indicates that the time averaged loss on the airfoil in an unsteady environment is substantially higher 
than it would be for the same airfoil in a steady environment. In addition, the minimum and maximum 
losses in this experiment show good agreement with calculated results by utilizing transitional and 
fully turbulent boundary layer prediction models. The increase in time averaged loss can be related 
to reduced frequency (Schultz (1977), Sharma et al. (1988)) as shown in Figure 2.4.7. Quasi-steady 
models developed on the basis of these results yield excellent agreements with the steady and unsteady 
time averaged boundary layer data from Hodson (1983) and Pfeil & Herbst (1979) as indicated in 
Figure 2.4.8. It should be emphasized that a relatively simple model can yield good predictions for the 
effect of unsteadiness on airfoil boundary layers once the primary physical phenomenon has been 
identified through review of basic experimental data. 



LOCATION Of UfSTNEAM NOTON NLAOC <KNCEMT WO* OF THE 2ND STATON > 

Figure 24.6 Mid— Span Losses on the 2nd Stator as Influenced by the 
Unsteadiness Generated by the Upstream Rotor Airfoil Wake. 
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Figure 24.7 Additional Time— Averaged Loss Generated Due to Unsteadiness 
Induced by Upstream Wakes can be Related to Reduced Frequency; 
Schultz (1977), Sharma et aL (1988, 90). 
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Figure 24.8 Quasi-Steady Boundary Layer Code Gives Good Estimate of 
Extra Losses Generated Due to the Effect of Unsteadiness. 

The stable of 2D CFD codes for turbomachinery design application were completed with the 
development of 2D Navier- Stokes code (Davis et al. (1986)). This code permits realistic simulation 
of airfoil performance at off— design conditions where foil separation effects become important. An 
example of the application of this code to estimate loading distribution on an airfoil is shown in Figure 
2.4.9 where the airfoil has a separation bubble on the airfoil pressure side. The Navier-Stokes code 
yields better agreement with experimental data than a 2D Euler flow solver because it simulates the 
effect of the separation bubble explicitly. It should be pointed out, however, that a good transition 
model is essential if good estimates of airfoil loadings are to be realized. 



Figure 2 4.9 2-D Navier- Stokes Code With Improved TYansition Model Shows 

Separation Bubble on Airfoil Pressure Side and Predicts its Effect on 
Airfoil Loading. 
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2.43 Second Generation of CFD codes for Tbrbomachinery 


Three-dimensional Euler codes (Ni et al. (1989), Denton (1975)) and Navier— Stokes codes 
(Hah (1983), Rhie (1986)) represent the second generation of CFD codes utilized in the 
turbomachinery design process. Application of these codes represents the state— of— the— art in 
turbine design and analysis systems. Availability of 3D CFD codes to analyze flow through turbine 
airfoil rows has greatly enhanced the capabilities of the design engineer. In addition to providing the 
ability to compute flow through complex 3D geometries, these codes also permit calculation of 
secondary flows in turbine airfoil rows and thus provide a vehicle to control these flows; this can result 
in significant improvement in performance and durability of turbines. 


The secondary flow prediction capability of the 3D CFD codes is demonstrated in Figure 2.4.10 
where experimental data from Langston’s Cascade (1977) are compared to theoretical predictions by 
utilizing both a Navier— Stokes code (Rhie (1986)) and an Euler code (Ni et al. (1989)). Both codes, 
run with measured inlet flow profiles, show good agreement with data in terms of airfoil loadings and 
surface streamline patterns, indicating the effects of secondary flows. The Navier— Stokes code was 
expected to yield good predictions for the three-dimensional flows since it solves viscous flow 
equations. Realistic predictions obtained tty using an invisdd Euler code were, however, surprising 
because secondary flow behavior was expected to be induced by viscous effects not modeled by the 
code. Further investigations suggested that good agreement between data and predictions due to 3D 
codes was mainly influenced by the use of measured upstream flow profiles; the predicted magnitude 
of secondary flows changed dramatically as the upstream profile changed in both codes. Specification 
of upstream boundary condition thus becomes an important variable in the turbomachinery design 
and analysis process. 



Figure Z4.1 0(a ) Both 3-D Euler and 3-D Navier-Stokes Code Yield Good 
Agreement with Langston ’s Cascade Data. 


U DATA (SUCTION SUNT ACC TI 



Figure 14.10(b) Both 3-D Euler and 3-D Navier-Stokes Code Show Good Agreement 
with Langston’s Cascade Data for Suction Surface Streamlines and 
Loading Distributions on the Airfoil Surfaces. 

Current turbines have closely spaced rows of airfoils. Tb analyze the flow through an airfoil row 
inside a multistage turbine, one must specify the boundary condition for that row which is difficult to 
obtain. Tb overcome this difficulty, multistage Euler codes have been developed (Ni et al. (1990) UJ 
these codes the flow downstream of each airfoil row is averaged at the interface plane (Figure 2.4.1 ) 
before information is transferred from an upstream to the downstream row during each solution 
iteration. All interstage boundary conditions are thus automatically provided by the code. Theoretical 
predictions from this code are compared against data measured in a two- stage turbine. Good 
agreement between data and predictions confirms the validity of the model. The measured airfoil 
loadings for the same experiment are compared to predictions in Figure 2.4.12 . Once again, excellent 
agreement is shown between the measured data and theoretical predictions. These multistage codes 
are routinely used in the design and analysis of the aircraft and SSME turbines. 
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Figure 24.11 Schematics of a Two— Stage TUrbine Showing the Strategy Used in 

Computing 3-D Flows by Using Multi-Stage Euler Code. Flowfield 
Downstream of Each Airfoil Row is Averaged and Calculations are 
Conducted in Their Frame of Reference. 
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Figure 24.12 Three— Dimensional Multi— Stage Euler Code Gives Good Estimates 

of Time -Averaged Loading for Airfoils in Unsteady Environment. 

The multistage calculation of Ni et al. (1989) described above facilitates steady state interaction 
between airfoil rows and provides good estimates of airfoil loadings, flow distributions and gap 
averaged total pressure and total temperature fields. 
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2.4.4 Future Direction 


The 3D CFD codes described in the previous subsection provide necessary tools for designers 
to optimize turbomachines under “mean” (steady) flow conditions. However, as pointed out earlier, 
the actual flow fields in turbomachines are highly complex and unsteady. The airfoils undergo 
transient and periodic aerodynamic loads which affect life and performance. In order to further 
improve the durability and efficiency of the machine, more advanced flow prediction codes need to 
be developed. 

TVvo specific aspects of the turbomachinery flow fields need further development: improved 
predictions of the viscous flow effects and a realistic account of the effects of unsteadiness. 


The viscous flows in turbomachines are dominated by complex boundary layers having laminar, 
transitional and turbulent flow regimes. There are currently no methods available that can predict the 
breakdown of laminar flows into transitional, and subsequently turbulent flows without utilizing 
empiricism. Development of models that can account for mixing induced by turbulent flows from first 
principles, although highly desirable, is beyond the scope of the present program. Even if such a 
method could be developed, it would predict a flow structure similar to that shown in Figure 2.4.13(a) 
which shows a flow field downstream of the rotor in the absence of upstream distortions. The real flow 
in the rotor is influenced by wakes and vortices from the upstream stator airfoils (Figure 2.4.13(b)). 
The overall flow structure in Figure 2.4.13(b) is totally different from that in Figure 2.4.13(a), 
j pHiVating that the effect of unsteadiness is more pronounced than the effect of viscosity. Further 
improvements in the turbomachinery flow field analysis are thus likely to come from unsteady flow 
analyses. 
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Figure 24.13(a) Rotor Exit Flow Total Pressure Loss Contours in the Absence of Upstream 

Wakes Show Distinct Organized Structures Similar to the Ones Measured 
in Cascades A 3-D Steady Navier- Stokes Code With Best TUrbulence 
Model Can Only Reproduce These Flow Structures. 
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22 Vanes 



Figure 14.13(b) Rotor Exit Flow Total Pressure Contours As They Are Influenced 
By Wakes From Upstream Stator Row Airfoil Show Large Pulsation 
in the Organized Flow Structure. To Predict These Flows, Viscous 
Unsteady 3—D Code is Required. 

Progress in computing three-dimensional unsteady flows has been reported by Rai (1987,89) 
where an innovative numerical scheme and a supercomputer (CRAY— XMP) were utilized to 
calculate three-dimensional viscous unsteady flows through a large scale model turbine. Some of his 
results are shown in Figure 2.4.14. Although the computational resources required for this exercise 
are large, Rai demonstrated that unsteady turbomachinery flow calculations are within the grasp of 
turbomachinery designers. Even though the geometry was simplified by assuming the same number 
of rotor and stator airfoil counts, the computation took 100 CPU hours on the CRAY-XMP 
computer. It is estimated, however, that a realistic case would require 1250 hours of CPU on the 
CRAY —XMP computer. Such requirements are not practical during a normal turbomachinery design 
process. In order to reduce these requirements, “smart” use of CFD codes must be adopted. Using 
the same approach as in the previous phases of the design system development, physical models, 
representing salient features of the flow field physics, must be developed from experiments. These 
models can then be incorporated into the unsteady, 3D viscous flow solvers, thus achieving improved 
accuracy with much reduced computational costs. 
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Figure 2 4.14(a ) Rai’s Unsteady Navier- Stokes Code Shows Good Agreement With 
Dring’s LSRR Stator and Rotor Airfoils Time-Averaged Loading 
Data. Similar Predictions Have Been Obtained by Using 3-D 
Multi-Stage Euler Code (Ni et aL (1990)). 
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Figure 14.14(b) Comparison of Pressure Amplitude Data on the Stator and Rotor 
Airfoils at Mid— Span Between Rcu’s 3—D Unsteady Navier— Stokes 
Predictions and Dring’s Experimental Data. An Improved Agreement 
Between Data and Predictions Indicated by Rai When He Used 3 
Stator and 4 Rotor Airfoils in His 2—D Unsteady Navier— Stokes 
Simulation. In Experiment, Rotor and Stator Consist of 28 and 22 
Airfoils Respectively. 

Work conducted in the present program has been a step towards developing the third generation 
of CFD code application in the turbomachinery design process. The development and verification of 
a 3D unsteady viscous flow solver in a turbomachinery environment is essential before information 
from these codes can be utilized to improve durability, structural integrity and performance of rocket 
engine turbopumps. Theoretical treatment including a system of equations, boundary conditions and 
turbulence modeling used in the CFD code development is discussed in Section 3. Verification of the 
code against benchmark quality data is discussed in Section 4. The verified code is applied to two 
configurations to provide accurate estimates of flows in radial and axial machines in Section S. 
Conclusions and recommendations for future work are discussed in Section 6. 
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3. THEORETICAL ANALYSIS 


The CFD code used in the present program was originally developed by Rai (1987, 89) to 
ci mniate unsteady two- and three-dimensional viscous flows through a turbine stage. Detailed 
discussions of the governing equations, integration procedure, turbulence model and boundary 
conditions are given by Rai and Dorney et al. (1992) and are reproduced here for the sake of 
completeness of documentation. Modifications to turbulence models and boundary conditions 
developed under the present program are also discussed. 

3.1 Governing Equations 

The code is based on numerical solution of the time dependent, three-dimensional, Reynolds 
averaged Navier— Stokes equations. These equations can be written in Cartesian coordinates as. 

Q, + (F, + F,\ + (G, + G,) y + (H; + H,) t = 0 (3.1.1) 


where 



QU 


0 

QU 2 + P 



guv 

- 


QU W 



(e, + P)u 



QV 


0 

QUV 


V 

QV 2 + P 

G, = - 


QV W 


** 

(e, + P)v _ 


. ^ . 

m 



QW 


0 

QUW 



QVW 

H,= - 


QW 2 + P 



( e, + F)w 




(3.1.2) 


(3.1.3) 


(3.1.4) 


(3.1.5) 


31 



and 


T, = 2pu, + X(u x + V M + w x ) 

** - P( u y + v ») 

r m - A«(«* + *x) 

W 

t„ -= + X(u x + V, + *V.) 

T„ “ M v « + **-,) 

T„ - T, 

T 9* T J« 

T„ - 2pw, + A(u x + V, + IV.) 
r fa * KT„ + vr* + HT. + ypP r ~'e x 
T v - MTj. + VTjy + MTj, + YflPr 'e, 
T fc *= ut„ + vt v + WT a + yfip-'e, 

e m e(y - i) 

*.* + «fe!±£±i3 


(3.1.6) 


For the present application, the second coefficient of viscosity is calculated using Stokes* 
hypothesis, X. = -2/3p. The equations of motion are completed by the perfect gas law which takes the 
form 


P = qRT (3.1.7) 

It is useful to non— dimensionalize the equations of motion so that certain parameters, such as 
the Reynolds and Mach numbers, can be varied independently. The non-dimensional variables 
chosen in this investigation are 
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In addition, for the analysis of arbitrary geometries the equations of motion can be generalized 
by vsi ng body— fitted coordinates. Using the independent variable transformations 

T = t 

£ - Z(x,y,z,t) (3.1.9) 

n - n{x,y,z,t) 

£ = £(*,**, 0 

the body— fitted Cartesian coordinates in the physical domain become uniform coordinates in the 
computational domain. The transformation allows easier implementation of boundary conditions 
since the geometry surface lies along the boundaries of the computational domain. Upon applying the 
transformations and non— dimensionalizing, the three-dimensional Navier— Stokes equations can 
be written as (where the superscripts have been omitted for clarity) 


Q r + (F, + Re _1 F„) { + (G, + Re"'G,), + (//, + R - 0 


(3.1.10) 


where 

F, (G,£) = /"'(££ + £F, + £,G, + £^) ( 3 . 1 . 11 ) 

G, {Q,ri) = J-\r\Q + nfi + V,G i + V#,) 

h,(Q,S) = J-\W + U 7 . + £A + wo 

The viscous fluxes are simplified by incorporating the thin layer assumption (Baldwin and Lomax 
(1978)). The thin layer assumption states that for high Reynolds number flows, the diffusion terms 
normal to a solid surface will be greater than those parallel to the surface. In the current study, viscous 
terms are retained in the direction normal to the hub/shroud surfaces ( 5 — direction) and in the 
dir ection normal to the blade surface (rj— direction). Thus, the non— dimensionalized and 
transformed equations now become: 

Qr + (F.)t + (G. + Re_1 G.), + (Hi + “ 0 ( 3 . 1 . 12 ) 


where 
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(3.1.13) 


(3.1.14) 
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The vector H, is obtained by replacing rj with 5 in Eqs. (3.1.13) and (3.1.15). If all the viscous 
terms are neglected, then the equations become the invisdd (Eider) equations of motion. The 
Jacobian of the transformation and the other metric quantities are given by Rai ( 1987, 89) and Pulliam 
and Steger (1985). 


J~ l ■= xpp t + xyp, + xyp t - X#p, - Xfft - xyp t 


(3.1.15) 
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Axy^-yfy) 

k = Ay& ~ w) 
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(3.1.16) 


The metric derivatives are evaluated using three point central differences in the interior of the 
computational domain and three point backward differences on the boundaries. 

3.2 Integration Procedure 

The governing equations of motion are integrated in time using the Approximate Factorization 
(AF) implicit technique developed by Beam and Warming (1977). Applying the AF technique for 
three-dimensional problems is accomplished by solving three one-dimensional operators, each 
requiring the inversion of a block tridiagonal matrix system with 5x5 blocks. Newtown iterations are 
applied within each global time step to increase stability and eliminate the linearization errors caused 
by the factorization process, lb apply Newton’s method, one starts with an initial guess for the solution 
and iterates according to: 


Qr + » = Qf - 

f{Q) 

This method can be applied to the unsteady Navier— Stokes equations by setting 


(32.1) 


f(Q) - G, + (F,) t + (G. + Re-‘G,) f + (//, + Re-‘tf,) { 
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where A is a Jacobian matrix, the factored, iterative, implicit integration algorithm can be defined by 
(Rai (1987, 89)) 


[/ + f§(v* w + v^)J [' + + 
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where 



(3.2.6) 


(3.2.7) 


and A, V, and 8 represent forward forward, backward, and central difference operators. In equation 
(2.21), & is an approximation to g" +1 . The quantities F* G h H h G„ and H, are numerical fluxes which 
are consistent with the physical fluxes F„ G„ H it G„ and //,. If p = 0 then & * O', and when the 

equation is iterated to convergence Qf - Q" +i . As the left hand side of equation (3.2.5) is driven to 
zero, the linearization and factorization errors associated with the AF technique are also driven to 
zero. If only one iteration is used then the integration scheme reverts to the conventional AF type 
scheme (Beam and Warming (1977)). Typically, unsteady calculations require two to t^ee iterations 
per global time step to reduce the residual of the density by three orders of magnitude (Rai (1987, 89), 
Domey et al. (1992)). 
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The inviscid numerical fluxes F it G h and H, are discretized using Roe’s scheme (Roe (1981)). 
The numerical fluxes are then evaluated from a family of high accuracy representations of the fluxes 
developed in Chakravarthy and Osher (1985). For example. 
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(32.8) 


AF + , 


-i/VJt 


where (F ( ) 1+1/Wk is the first order accurate upwind flux given by 

(F i),+i/i^> * (Fi)i+i^jk] ” 2 ^^ ^ i+uv*) (32.9) 

and the additional terms in Eq. (32.8) are used to increase the order of accuracy. The third order 
accurate upwind biased difference scheme is used for interior grid points and either first or second 
order accurate upwind differencing is used at boundary points. 

The flux differences (AF 1 ) in Eqs. (3.2.8) and (3.2.9) are calculated using Roe’s scheme and are 
given by 


AF* j+1/w = A 4 i+1/a * x (& +1 * - Qm) (32.10) 

The flow variables needed to determine A 4 between grid points (i + 1/2, j, k) are calculated using 
Roe’s averaging formulae: 
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where h t is the total enthalpy and is defined as 


h, = 


*,+/> 
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(32.12) 
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The Jacobian matrix A * can be written as: 

A ± = T ( ~ l AfT ( 

where T t and Tf are eigenvectors and which can be expressed as (Pulliam and Steger (1985)): 
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The eigenvectors for the B and C matrices can be obtained by replacing | in Eqs. (3.2.14) and (3.2.15) 
with Y| and 5 . The matrix containing the eigenvalues of the Jacobian matrix is given by 




± 





(32.17) 


where 


- £, + Im + Z,v + £ t w 

X 4 « A, + xa (32.18) 

X s * J, — xa 


The superscripts ‘±’ in the previous equations refer to the contributions from the downstream and 
upstream running characteristic waves, lb prevent expansion shocks, the eigenvalues in Eq. (3.2.18) 
are replaced by a nonvanishing, continuously differentiable approximation which can be written as 
(Yee et al. (1985)): 


where n = 1,5 



|A.| < » 
|J„| * 6 


(3.2.19) 


The flow variables needed to determine the viscous fluxes, G. and H„ are evaluated using 
standard central differences. For example. 


(^»)v + l/2Jk = /(Gy + l/2Jk.(G»)v+lA») 

Qii+i/vt “ 2 ^** + Cy+w) (3.2.20) 

(Qn)ij + l/2Jt = Qij + ljc Qijjt 


The corresponding viscous flux Jacobian, M, ran be written as (Pulliam and Steger (1985)) 
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where 


m 21 = - 5,a,(u/e) - Sjd,(v/e) - s 3 d,(w/g) 
m s , — - Sjd,(«/e) - S t d 4 (v/g) - Sjd,(w/g) 
m 41 - - 5,a,(«/e) - Sfjy/e) - S t d,(wjg) 
m 5l - - S„d,[- We 2 ) + (« 2 + y2 + "Vel 

- s,a,(uVe) - Sjj?/Q) - s t d,(w 2 / e ) 
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The viscous flux Jacobian, N can be obtained by replacing tj in Eqs. (3.2.21) and (3.2.22) with 5 . The 
equations of motion and solution procedure used in the two-dimensional computational procedure 
are a direct subset of the equations developed above, except that the inviscid fluxes are calculated 
using Osher’s (1981) approximate Riemann solver. 

33 Tbrbulence/TVansition Model 

The original version of the CFD code used in the present program contained a simple turbulence 
model which facilitated reasonable predictions in an unsteady flow environment. Limitations of this 
model were that it assumed fully turbulent flows on airfoil and endwall surface whereas in reality flow 
is transitional on airfoil surfaces especially in available benchmark quality experiments against which 
this code was planned to be verified. Appropriate modifications were, therefore, made to the 
turbulence model to allow transitional flow calculations. Additional modifications were also made to 
account for free stream turbulence level, surface roughness and extra rates of strain to allow realistic 
predictions in the complex flow environment of the rocket turbopump. All of these modifications were 
maHf. in the context of the base turbulence model, these modifications are, however, of generic nature 
and would still be applicable if a higher order turbulence model was implemented in the code. The 
scope of the present program limited implementation of the higher order turbulence model in the 
code which is needed to obtain further reliable predictions from the code. 
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An eddy viscosity formulation is used in the code to model the effect of turbulence. Effective 
viscosity and effective turbulent conductivity are defined as: 

Effective viscosity = fi = fi L + n T (3.3.1) 

Effective conductivity = (3-3.2) 

I'P **1, n r 

Ptl and Ptt are the laminar and turbulent Prandtl numbers respectively. 

33.1 Baldwin — Lomax lUrbulence Model 

Ibrbulence model developed by Baldwin- Lomax (1978) was used in the original CFD code used 
in the present program. 

In the Baldwin-Lomax (B-L) model turbulent eddy viscosity \ij is described by 





s £ s, 
s > s, 


(3.3.3) 


where s is the distance normal to the solid surface and Scronover is the smallest value at which 
fi T . — n T . In the inner region, the eddy viscosity is calculated using the Prandtl— Van Driest 

formulation 
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and the magnitude of the vorticity, | <d | can be written as: 
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In the outer region the eddy viscosity is calculated using 

Ht = KC^qF bJuF nufe) 

where K is the Clauser constant, Cq> is an additional constant and Fwake is described by 

«= iam(s tu J 7 mm C ¥lt s ma 4jf/F mn 


(33.8) 

(33.9) 


The term F ma * is the maximum value of F(s) along a given computational grid line normal to the 
surface and 


F(s) *= s\w\D 

The Klebanoff intermittency factor, Fncb(s) is given by: 

F»Js) - (1 + 5.5(($C ja *)/s n «) # ) 


(3.3.10) 


(3.3.11) 


and qdif is the difference between the maximum and minimum velocity in the profile. The vorticity and 
velocity are calculated in the appropriate reference frame (i.e., in the absolute reference frame for 
stationary surfaces and in the relative frame for moving surfaces). The constants used in the B L 
turbulence model are: 


A* - 26 
Cm, = 03 
k = (14 


C„ » 1.6 
C* = a 25 
K = .0168 


(33.12) 


The Baldwin-Lomax turbulence model is based upon two-dimensional boundary layer data 
and as such, is not well suited for corner flows such as those at the blade/endwall juncture. 

Originally, the treatment used to implement this turbulence model in the corner regions (Rai 
(1989)) was the technique proposed by Hung and Buning (1984). In this technique, the turbulence 
model is computed separately for each endwall and the blade surface. The mixing length m the corner 
region is computed depending on the computational indices of a given node. For instance, consider 
the case when the J= constant computational lines run normal to the blade and the K-constant lines 
run normal to the endwall. For any computational node whose J-wise index is less than ils K-vnse 
index, the normal distance is defined as the distance from the blade surface to the grid point and the 
parallel distance is defined as the distance from the endwall to the grid point. The mixing length for 
the inner region of the boundary layer is then calculated as 

/ = 2sn/(s + n + /(r 2 + n 2 )) (33.13) 

where s is the parallel distance and n is the normal distance. The eddy viscosity is then based on die 
flow variables along a computational grid line from the airfoil surface to the grid point under 
consideration. Likewise, for any computational node whose J-wise index is greater than its K wise 
index, the parallel distance is measured from the blade surface to the grid pomt and the normal 
distance is measured from the endwall to the grid point. The eddy viscosity is then based on the flow 
variables along a computational grid line from endwall to the grid point. TWo significant problems 
arise from this particular three-dimensional implementation of the Baldwm-Lomax turbulence 
model. First, the eddy viscosity distribution in the corner regions is discontinuous across the J 
computational lines and can cause large gradients to occur in the velocity field. Secondly, this 
particular blending is dependent upon the computational grid density and stretching in both 
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directions. It was found, based upon numerical simulations, that flow solutions in the blade/endwall 
region were extremely sensitive to changes in the computational grid structure. 

A blending function was then used to smoothly vary the eddy viscosity distribution between the 
blade and endwall. Separate eddy viscosity distributions are computed for the blade and endwall 
surfaces along the computational lines which run normal to each surface, respectively. The eddy 
viscosity in the corner flow regions is then computed based upon the following blending function 
according to the work of Vatsa and Wedan (1988) 


Mr “ 


tlMr k + tl M Tm . 

il + li. 


(33.14) 


where lb is the distance from the blade surface to a given node, l«w is the distance from the endwall 
surface to the node, and fx-rb and are the eddy viscosities computed from the separate blade and 

endwall flows, respectively. This type of blending creates a smooth eddy viscosity distribution in the 
comer regions. 


333 Turbulence Model for Surface Roughness 

The need to model surface roughness effects was illustrated by tests conducted at NASA MSFC 
for the SSME turbine (Boynton et al. (1992)). Significant improvement in the performance of the 
turbine was achieved by polishing the hardware currently used in the SSME. Surface roughness erodes 
the effect of viscous damping near the wall causing an increase in mixing length in the inner part of 
the boundary layer. A model proposed by van Driest (1956) can be used to account for this 
phenomena. The modified damping function due to surface roughness is given as: 

Dm = D + Dr (3.3.15) 

Dm = damping term for modifying turbulent viscosity 

D = damping term for smooth surface (equation 3 33) 

Dr = -60 y + /r + A + 


y + is defined in equation (33.7) 

A + is defined in equation (33.12) 

r + can be obtained by using surface roughness height instead of ‘s’ in 
equation (3.3.7) 

Predictive capabilities of this model were demonstrated in a publication by Blair (1992) through 
comparison against data obtained on a smooth and rough airfoil as a part of the NASA MSFC contract 
# NAS8— 37351 as indicated in Figure 33.1. 
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D Smooth, K * 7.6 micron* 

• N**r-Smooth, K - SI micron* 
A Rough, K - 660 micron* 


STANS Prediction* 
Roughn*** valo* indic*t*d 
In micron* 



S/Bx 


Figure 3.3. 1 Comparison of Midspan Rotor Airfoil Heat Thansfer Distributions 

Obtained at Re = 5.8 x 10 s and fa = 40 * with STAN5 Predictions 
for Various Wall Roughness Values. 

333 Turbulence Model for Extra Rates of Strain 

Extra rates of strain such as streamline curvature, system rotation and streamtube contraction 
ratio have been known to have significant effects on the structure of turbulence as discussed by 
Bradshaw (1973), Johnston (1970) and Sharma and Graziani (1982). A solution of Reynolds stress 
transport equations is actually needed to resolve these effects properly, but the development of 
closure models for these equations is still in infancy. Simpler models are, therefore, needed to capture 
the effects of these flow features for engineering application. Modifications to turbulence model 
developed in the present program are based on the work of Sharma and Graziani (1982), who utilized 
Tbwsend’s hypothesis (1956) which states that the turbulent kinetic energy and Reynolds shear stress 
are related by a constant in the inner part of the turbulent boundary layer. One can, therefore, capture 
these effects of extra rates of strain by interrogating the transport equation for turbulent kinetic energy 
and the total Reynolds shear stress. Modified turbulent viscosity can be fairly accurately estimated as 
the ratio of the generation terms for the Reynolds shear stress equation to the production terms for 
the turbulent kinetic energy transport equation. No empirical constant is required in this approach 
and it yields fairly good results for flows in the presence of extra rates of strain. 

For incompressible flows production terms for the three normal stresses and the turbulent 
kin etic energy are given below in Cartesian coordinate with rotation in the axial, tangential and 
normal direction: 


+ (33.16) 

(3.3.17) 
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+2fl, )“ w (^" 2 ^) (3 - 318) 

Adding the above three equations yields production terms for the turbulent kinetic energy 


(*-?+?+? > : 


- u * W _ rrdV_-idW 

“a* v dy w dz 


— (¥♦«)— (f + f) < 3JI9 > 


It is apparent from the above four equations that whereas each individual normal stress 
component is affected by rotation (Qx, y, z), total turbulent kinetic energy is unchanged indicating that 
turbulence models based on turbulent kinetic energy transport equation would fail to capture the 
effect of rotation explicitly. Since rotation terms change the distribution of energy in various 
components, they could either enhance or reduce mixing and, therefore, losses and heat loads in 
turbomachines. 

Generation terms for the diagonal components of Reynolds shear stress are given below to show 
a more complex effect of rotation in turbulent flows, effects of streamtube contraction and streamline 
curvatures are also present in these terms but these are not apparent due to the use of the coordinate 
system: 

G(0P) - + + £) + p- + 2^) + P KW + “*) + m (%- **) 

(3.3.20) 

Cm = fl*(f + f) + «^f + f) + P- - Zf-g + 2Qy) + + 2ftc) 

(3321) 

^ + f ) + P " + **) + + H + ' “*) 

(33.22) 

A generation term for the magnitude of the diagonal stress (defined as t — /op 3 + tw 2 + uw 2 ) 
can be deduced from the above three equations by multiplying equations (33.20), (3321) and (33.22) 
by and pw/t, respectively, as shown below: 


G( t) - TERM (D + TERM @ + TERM Q> 


(3323) 


where 


term ® . + 20.) ♦ + 20,) + ♦ 20.) 

n* O - ?(f ♦ f ) ♦' ¥(g ♦ f) ♦¥(£' - f ) 
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TERM ( 3 > in the above equation represents a generation of Reynolds shear stress due to the 
effects of rotation and curvature. This turbulence generation mechanism is not present in the 
turbulent kinetic energy production terms as indicated by the absence of rotation terms in equation 
(33.19). 

TERM CD represents the generation of Reynolds shear stress due to the effect of streamtube 
convergence and divergence. 

TERM (D represents the direct generation of Reynolds shear stress. It was pointed out by 
Bradshaw (1968) that the generation of Reynolds shear stress is also affected by the pressure-velocity 
correlation terms present in the complete Reynolds shear transport equation. In the absence of extra 
rates of strain for equilibrium turbulent boundary layer flows, generation of turbulence given by 
equation (3.3.19) and by TERM ® modified to account for pressure-velocity correlation have to be 
the same, therefore, direct generation of Reynolds shear stress can be written as: 

TERM (D + Extra terms from pressure - velocity correlation = jrP(K) (33.24) 

Generation of the diagonal component of Reynolds shear stress can now be written as: 

G(t) = ±P(K) + (TERM © + TERM (D) 

(33.25) 

= L Pf J i + K (re*A/©+ TERM (D) 1 

T P(K) J 

This equation shows that production terms in the turbulent kinetic energy equation or other 
turbulence mixing equations must be modified to properly account for the effects of extra rates of 
strain. These modifications have invariably been conducted by using ad hoc expressions, present 
analysis indicates that these can be exactly deduced from turbulence transport equations. 

In equilibrium turbulent boundary layer flows, production and dissipation terms in the turbulent 
kinetic energy transport equation balance each other. Dissipation and generation terms in the 
Reynolds shear stress transport equation must, therefore, also balance each other in equilibrium 
boundary layer flows. Using Bradshaw’s (1968) model for dissipation terms, the behavior of Reynolds 
shear stress transport equation can be expressed as: 

G( t) = (Dissipation) 

(33.26) 

T W 2 

~ K L 0 

where 

L 0 = dissipation length 
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By substituting equation (3.3.26) in equation (3.3.25) results in: 


(TERM @ + TERM ® 


)[i+f 




- 

I/O 


+ TERM® 

m 


In the absence of TERM Q) and TERM Q> , it can be shown that 


I- - "t(i£ + H) 


(3.327) 


(3328) 


(3329) 


Hr in this equation can, therefore, be modified to account for TERM Q) and TERM (D as: 

Hi „ - Hr (l + ¥ *0.* + ¥ *») <3 ' 3 ' M) 


where 


H Tm « modified turbulent viscosity 

ft T - turbulent viscosity from Baldwin — Lomax model discussed in Section 33.1 

and TERM (5) 

RiQ,R « Richardson number for rotation & curvature = — — 

TERM ® 

/b'jo = Richardson number for three — dimensionality = — — 

Impact of these modifications can be illustrated in boundary layer flows for two-dimensional 
rotating ducts and at the line of symmetry in cascades. 

In two-dimensional boundary layers experimental data (Klebanoff (1954)) indicates that 
£ - u? - K,^ - 0l 4K and iv* - (I6K 

In two-dimensional rotating ducts, equation (33.26) simplifies to: 


f*T m * f*T 


( 1 + f f) 


which is almost exactly the same as deduced by Johnston (1971) from experimental data. 

For boundary layer flows at the line of symmetry equation (3326) simplifies to: 

( dU'X 

‘-¥ f ) 

which was used by Sharma & Graziani (1982) to yield good estimates of heat loads on the suction side 
of a turbine cascade at the mid-span as indicated in Figure 332. 
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Modification to turbulent viscosity suggested by equation (3.3.26) are generic and can be 
implemented in one— or two— equation turbulence models almost exactly in the form presented by 
the equation. 

(a) (b) 



Figure 3.3.2 TUrbulence Model Modified to Account for Extra Rates of Strain 
(Sharma & Granani (1982)) Yields Good Estimate of Midspan 
Stanton Number on Langston’s TUrbine Airfoil 

(a) Nominal Inlet Boundary Layer to the Cascade. 

(b) Thin Inlet Boundary Layer to the Cascade (Refer to Figure 2.3.2 
for the Magnitudes of Flow Convergence at the Cascade Airfoil 

Midspan). 
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3.3.4 Transition Model 


Although all boundary layer flows start in a laminar mode, they become turbulent by passing 
through a tr ansi tional region, the regions of laminar and transitional flows are only significant in 
turbomachinery operating under a Reynolds number of about one million. For rocket applications 
Reynolds numbers are invariably larger than ten million and these applications transition plays a small 
role on the overall development of the viscous flows. A knowledge of the transitional region was, 
however, essential during the code verification and application part of this program since most of the 
available data used in these tasks were acquired at moderate Reynolds numbers where laminar and 
transitional regions were significant. Transition models were, therefore, implemented in the CFD 
code to accomplish the defined tasks. 

In three-dimensional flow situations, the transition from the laminar to the turbulent state was 
assumed to be instantan eous. The transition point was specified along grid lines to simulate the 
physics of the flowfield. 

For two-dimensional and unsteady flow computations the onset of transition was determined 
by using a correlation based on the work of Mayle (1991) and Hourmuziadis (1990). According to this 
correlation, onset of transition occurs when the local boundary layer momentum loss thickness 
Reynolds number exceeds a critical Reynolds number. 

R0c = 40//fu (33.24) 


where 

R6c = Critical Reynolds number, 
and 

TU = Local turbulence level which needs to be specified as input to the computer 
code for single airfoil row or for the first row of the stage, it is automatically 
calculated for the following rows. 

Numerical computations are conducted by modifying the turbulent viscosity in the transitional 
region by using the approach suggested by Sharma (1987). 


UTM 

Ptm 

Ut 

F 


R0 


F*|ij 

Modified turbulent viscosity 

TUrbulent viscosity calculated in Section 33.1 

Intermittency factor having a value of zero in laminar regions 

and a value of 1 in fully turbulent region 

1 - exp ( - (R0 2 - 3 - ROc^/ROc 2 - 6 *) 

Reynolds number based on momentum loss thickness of the 
boundary layer. 


(3335) 
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3.4 Boundary Conditions 
3.4.1 Inlet Boundary Condition 

In the following expressions, u, v, w are velocities in the x— , y— , and z— directions, p is static 
pressure, q is density, c is speed of sound, s is entropy, and h is total enthalpy, ci is arctan(v/u) and 4* 
is arctan(w/u); boldface type indicates a new updated value, subscript 0 indicates an initial value, and 
subscripts of 1 and 2 indicate a current value at inlet and one point downstream of inlet, respectively. 

The following boundary conditions were available in the original version of the Rai code (Rai 
(1989)): 


Reimann invariant 1 = R t = u 0 + ^ ^ j— 

| (3.4.1) 

2 ( Jypl\ 

Reimann invariant 2 - R 2 - u 2 J 

| (3-4.2) 

«i = (R x + R 2 )/2 


Vj = v 0 

(3.43) 

Wj = H-0 



(3.4.4) 

Sj — s 0 

(3.4.5) 



* = w 



(3.4.6) 

„ QSi 2 

P l " y 



The following modifications were made to, in effect, hold inlet total pressure, in addition to 
reducing the stiffness of the boundary condition: 

First define the variable h’ as follows: 


h' - + u 1 = 2h - (y 2 + w 2 ) 

y - 1 

(3.4.7) 

Assuming hi = ho, that is, inlet total enthalpy is held at its initial value, 

h'l = 2/l 0 — (V! J + H>, 2 ) 

(3.4.8) 

It can be shown that (3.4.8), together with the assumption 

2c, 2c, n 

y-l “ 2 y — 1 **’ 

(3.4.9) 

yields: 

(3.4.9a) 
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(3.4.10) 


This is a quadratic equation for ci that can be solved with the quadratic formula. Then, 


u 


= R 2 + 


j£i_ 

Y- 1 


v, - «, tan(a 0 ) 
h-, *= «, tan(4> 0 ) 


5, =5 0 


Pi 



Pi 


££i! 

y 


(3.4.11) 

(3.4.12) 


3^2 Surface Boundary Conditions 

The surface boundary conditions in Rai code (Rai (1989)) were a first-order adiabatic wall 
condition, that is, in addition to the no— slip condition for velocity, zero normal pressure and density 
gradients were imposed. Options were added to specify either wall temperature or heat flux on the 
surfaces, which modified the density condition. 

The specified wall temperature condition that was implemented was developed by Griffin 
(1990). It requires an additional input file containing adiabatic wall temperatures on the surfaces, and 
the wall temperature is prescribed as a percentage of the adiabatic wall temperature. The percentage 
is input in the RAUOB shell script. The adiabatic wall temperature file can be obtained by running 
RAI3DS with the adiabatic wall condition; the file will get written out in routine WALTMP, which will 
also read the file if the option to specify wall temperature is chosen. The UNICOS jobstream will need 
modification to assign these files to the file environment. 

The following is the derivation of the second— order density condition allowing specified heat 
flux. Here, Q is heat flux, k is the coefficient of thermal conductivity, R is the specific gas constant, 
T is temperature, p is static pressure, q is density, and y is distance from the wall. 

Assume 


Q = 



- k 


d (*l -*($)* -/>(♦) 

dy “ R e 2 


Assuming a zero normal pressure gradient, this yields 

n ]v_(de\ ^ 4 > e H2' 

u Re 2 \4y) *y p 


where Q' =^r 


Now assume that q is some quadratic function of y, given by 

q ■= ay 2 + by + c 

Then, values of q at yl, y2, and y3 are given by the system 


ay 2 + fyl + C = 0, 


ay 2 + by 2 + c = ei 


ay 2 + fyj + c « es 


(3.4.13) 


(3.4.14) 

(3.4.15) 


(3.4.16) 


£ 
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Assuming yi is zero, that is, at the wall, the system reduces to 

, _ ( 02 - gi)y» 2 - (gj ~ Pi))': 1 

y#* - ytis 


Differentiating the quadratic equation for Q, 


^ = lay + b 
dy 


^ = b at the wall . 
dy 


From the previous expression for 


de 
dy * 


tea - eifr> 2 - (p3 - giW = 


Algebraically solving for Qi , 


e^ 2 - Q#2 2 _ qWwI Uiy£ (3.4.17) 

e ‘ = y'-yi p(y> 2 ~y2 2 ) 

The values of p and q at y 2 are used in the implementation of the condition, lb define the value 
of Q' (note the solver will expect the given heat flux value to be properly nondimensionalized and to 
contain the factors of R and k), divide the dimensional Q by the dimensional k, and use a value of 1.0 
for R to be consistent with the solver’s nondimensionalization scheme. This will result in a value with 
units of degrees over length. Divide the resultant value by the inlet total temperature and convert the 
re maining length unit to inches to obtain the value of Q’ in inches. This is what should be input to the 
RADOB shell script as discussed by Belford (Appendix A). 
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4. CODE VERIFICATION 


The predictive capability of the numerical solution procedure (code) can be verified by 
comparing the results from the code with: 

• Results from exact calculation methods for an unsteady multistage 
turbomachine 

• A complete set of unsteady experimental data for a multistage turbomachine 

• A series of experimental data set from benchmark quality experiments which 
simulate various aspects of multistage turbomachines. 

Exact calculation methods for an unsteady multistage turbomachine are not yet available. In 
addition, a complete set of experimental unsteady data for a realistic turbomachine is also not 
available. The verification of the code in the present program was, therefore, conducted by comparing 
its predictions with data from benchmark quality experiments which simulated pertinent aspects of 
multistage turbomachines. These comparisons, discussed below, dearly demonstrate that the code 
provides accurate estimates of loadings, losses and heat loads for airfoil rows both in a steady and in 
an unsteady flow environment. This indicates that the code will provide accurate and reliable 
estimates of flow fields in a multistage turbomachine and it can be used in the design process to 
improve the performance and durability of turbopumps used in the rocket propulsion system. 

4.1 Verification of the 2D Steady Aspect of the Code 

Predictive capabilities of the code in a two-dimensional steady flow environment were verified 
by comparing theoretical predictions with experimental data obtained at the mean section of three 
linear cascades and two airfoil rows in a stage environment. Although flow through the airfoil rows 
in a stage environment is unsteady in a strict sense, simple assumptions are made to treat the flow as 
steady to demonstrate that steady flow assumptions do yield fairly accurate estimates of the 
performance and heat load characteristics of airfoil rows. This exercise indicates that appropriate 
assumptions permit cost effective evaluation of airfoil geometries in a multistage turbomachinery; 
computational resources required to execute a steady flow simulation are almost two orders of 
magnitudes lower than those needed to execute an unsteady flow simulation. 

Representative grids used in the 2D code verification are shown in Figure 4.1.1. Nominally they 
contained 101*21 points in the inner grid and 50*31 points in the outer grid. An extensive evaluation 
of the effect of the grid points on the accuracy of loss and heat transfer predictions was conducted as 
part of work funded by the Naval Air Systems Command under NAVAIR Contract 
#N00140— 88— C— 0677 at the United Technologies Research Center by Domey, Davis and Edwards 
(1992). Since the base CFD code in this effort was the same as that used in the present program, results 
from the NAVAIR Contract are directly applicable to the present program. The focus in the present 
program is to achieve engineering accuracy of the solution using minimum computer resources 
without compromising the technical results. 
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Future 4.1.1 A RA13DC Stator Grid 


4.1.1 Kopper’s Cascade 

The first 2D test case termed as Kopper’s cascade was tested in the United Technologies 
Research Center (UTRC) Variable Density Supersonic Plane Cascade Wind TUnnel. Tins cascade had 
an aspect ratio of 237 and pitch to chord ratio of 0.8925. The air entered the cascade with an inlet angle 
of 313 degrees and a Mach number of 0.4187. The Reynolds number based on axial chord (1368 
inches) and exit velocity was 500,000 and it was tested with a pressure ratio (exit static 
pressure/upstream total pressure) of 0.625. At the operating condition this airfoil had a separation 
bubble on the airfoil pressure side which affected the airfoil loading distribution on both the pressure 
and the suction sides. This configuration was used in the code verification process to demonstrate that 
accurate simulation of the transition process is essential to capture relevant features of the flow field. 
Theoretical calculations were conducted both in a transitional and a fully turbulent mode. 

Theoretical predictions from the code are compared against experimental data for airfoil 
surface static pressure distributions (Figure 4.13). Results from transitional calculations exhibit 
excellent agreement with the experimental data while the fully turbulent calculations miss 
experimental behavior in the middle 40% of the airfoil. 
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Figure 4.1.2 Airfoil Pressure Distribution for Hopper’s Cascade 

The nature of the flow near the surface of the airfoil can be more dearly highlighted through 
review of the streaklines around the airfoil streaklines generated from the fully turbulent calculations 
and the transitional calculations are shown in Figure 4.13. Both calculations predicted separation 
near the pressure side of the airfoil as indicated in this figure; however, the separation bubble 
predicted by the transitional calculation is larger and better defined, yielding better agreement with 
the pressure distribution data on both the pressure and the suction sides of the airfoil as indicated in 
Figure 4.13. 



Figure 4.1.3 Predicted Streaklines Near Hopper’s AirfoiL 

(a) Fully Ihrbulent Calculation, 

(b) Ihmsilional Calculation. 


4.13 Hodson’s Cascade 

The second 2- D test cascade termed as Hodson’s Cascade was tested at the Whittle Laboratory 
at Ca mb ridge University in UJC. This linear cascade contained seven airfoils with an aspect ratio of 
3.0 and a pitch to chord ratio of 0.698. Tfest conditions involved ambient air entering the cascade at 
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approximately 58.89 ft/sec. The Reynolds number based on axial chord (32 inches) and exit velocity 
was 315,000. Additional details are given by Hodson (1983). 

This airfoil had a small separation bubble on the airfoil suction side downstream of the maximum 
velocity point. Comparisons of predictions with experimental data for this cascade were made to verify 
the boundary layer and the transition prediction capability of the code. Calculations were conducted 
both in a transitional and a fully turbulent mode. 

Measured air foil surface static pressure coefficients are compared to theoretical predictions in 
Figure 4.1.4. Although both fully turbulent and transitional calculations show good agreement with 
the data, transitional calculations show slightly better agreement with the data on the airfoil suction 
side. 



Figure 4.1.4 Airfoil Pressure Distribution for Hodson ’s Cascade 

Measured boundary layer parameters on the airfoil suction side are compared to theoretical 
predictions in Figure 4.1.5. Plots of boundary layer momentum loss thickness are presented in Figure 
4.1.5a. Results from transitional calculations exhibit excellent agreement with the data, showing 
improvement over fully turbulent calculations. Figure 4.1.5b shows a comparison of predicted and 
measured shape factor (displacement thickness/momentum loss thickness). Again results from 
tr ansi tional calculations agree more closely with the data, although the transitional results appear to 
overshoot the data in the trailing edge region. This overshoot is a result of the integration method used 
in calculating the integral parameters, the displacement and the momentum loss thicknesses. The 
definition of these parameters calls for integration along lines normal to the surface. The grid lines 
are normal to the surface near the airfoil but become quite skewed away from the surface at the trailing 
edge. A better method of integration would be to interpolate the results from the calculations onto 
a grid which has normal lines at the trailing edge before calculating the boundary layer integral 
parameters. Such a method would provide even better agreement with the experimental data than that 
indicated in Figure 4.1-5. 
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(a) (b) 

Figure 4.1.5 Boundary Layer Parameters for Hodson’s Cascade. 

(a) Momentum Thickness, 

(b) Shape Factor 


4.13 Dring’s Stator (Midspan) 

The third 2D test configuration termed as Dring’s Stator (midspan) was tested at UTRC as the 
first row of the 1% stage large scale rotating rig (LSRR). Although flow in this airfoil row is affected 
both by three-dimensionality (due of radial pressure gradient) and unsteadiness (induced by relative 
movement of the downstream rotor), these effects are not very pronounced and a two-dimensional 
simulation yields fairly reasonable estimates of flow through the midspan of the airfoil. The aspect 
ratio for this annular cascade is 1.0118 and pitch to chord ratio at midspan of 13. Ifest conditions 
involved ambient air at approximately 75 ft/sec. The Reynolds number based on axial chord (6 inches) 
and exit velocity is 612,000. Additional details are given by Blair et al. (1988). 

Detailed heat transfer data are available at the mean section of this airfoil. Comparisons of 
predictions with experimental data for this configuration were made to verify the heat load prediction 
capability of the code. Calculations were conducted both in a transitional and a fully turbulent mode. 

Measured airfoil surface static pressure coefficients are compared to theoretical predictions in 
Figure 4.1.6 showing excellent agreement between data and predictions; both transitional and fully 
turbulent calculations yielded identical predicted values for static pressures. 
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Figure 4.1.6 2D Steady Version of the Code (RAI2DC) 

Developed to Verify the Code Against Basic Data 

Measured heat transfer coefficient (Stanton number) distributions along the midspan of the 
airfoil are compared theoretical predictions in Figure 4.1.7. Fully turbulent calculations are found to 
overestimate the heat transfer coefficients on the pressure side and over initial 60% of the suction 
side. The transitional calculations are, however, found to yield excellent agreement with the data over 
both sides of the airfoil. This figure illustrates need to model transitional nature of the boundary layer 
to accurately estimate heat loads on airfoil rows. 


data predictions 

AVAIUULE IMPROVED 

< TURBULENCE) (TRANSITIONAL) 



Figure 4.1.7 Transition Model in the Code (RAI2DC) Provides Improved 
Estimates ofDring’s Stator Airfoil Surface Stanton Number 

4.1.4 Dring’s Rotor (Midspan) 

The fourth 2D test configuration termed as Dring’s rotor (midspan) was tested at the UTRC as 
the rotor in the IV 2 stage LSRR. Flowfield on this airfoil is strongly affected by unsteadiness due to 
rotation of the rotor relative to the adjacent stators. Measured time-averaged experimental heat 
transfer and static pressure data on this airfoil, however, indicated that both of these parameters were 
unaffected when the rotor was placed at two distinct distances (15% and 65% axial chord) downstream 
of the first stator. This result showed that it was possible to simulate the flow through this airfoil by 
assuming a steady flow assumption. Only fully turbulent calculations were conducted because of the 
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presence of a relatively large inlet turbulence level generated by the upstream stator . The aspect ratio 
of this rotor is 0.9464 and pitch to chord ratio at midspan of 1.01. Tfest conditions involved almost 
ambient total pressure with an axial velocity of 75 ft/sec and a relative inlet angle of 40.6 degrees. The 
Reynolds n umb er based on axial chord (6 indies) and exit velodty is 525,000. Additional details are 
given by Blair et al. (1988). 

Measured airfoil surface static pressure coeffidents are found to be in good agreement with 
theoretical predictions as shown in Figure 4.1.8(a). Fully turbulent calculations modified to account 
for free stream turbulence level and surface curvature are found to yield better agreement with the 
airfoil surface Stanton number (heat transfer coeffident) data (Figure 4.1.8(b)) than the base 
turbulence model available in the original Rai code. This result indicates that steady flow simulations 
fan yield fairly reliable estimates of time— averaged loadings and heat loads on airfoil rows operating 
in an unstead y flow environment; appropriate modifications are, however, needed to the turbulence 
model to account for physical variables existing in the flowfield. 
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Figure 4.1.8(a) 2D Steady Cascade Version of the RAI Code ( RAI2DC) 
Developed to Verify the Code Against Basic Data 
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Figure 4.1.8(b) Improved TUrbulence Model Yields Better Agreement with Rotor 
Stanton Number Data Than Available Model in the RAI Code 
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4.1.5 Transonic Cascade 


The transonic flow prediction capability of the code was verified by comparing calculations to 
d at a ob tain ed for a Pratt & Whitney Transonic Cascade. This cascade, tested in the UTRC Variable 
Density Plane Cascade Wind TUnnel, had an aspect ratio of 4 and pitch to chord ratio of 0.647. The 
air entered the cascade at an inlet angle of 48.5 degrees and a Mach number of 0.48. ’Die Reynolds 
number based on axial chord (1.5 inches) and exit velocity was 760,000 and it was tested with a pressure 
ratio (exit static/upstream total pressure) of 0.528. Calculations were conducted both in a transitional 
and a fully turbulent mode. The shape of the airfoil and airfoil surface static pressure distributions are 
shown in Figure 4.1.9. Both transitional and fully turbulent calculations are shown to yield good 
agreement with the experimental data except in the trailing edge region. The inability of the code to 
predict pressures in the trailing edge region can be attributed to the thin layer assumptions used in 
the solution procedure. The calculations are, however, found to yield good estimates of lift on the 
airfoil in^i™tin£ that it would provide accurate predictions of flows around transonic airfoils. 


* (Miih Ml) „ . 

dim C«lcft«4i*a *ii»| irrlsUrt* aaMI 

(AIMS Ctlcalsli* t I'M u-kilattt m««I 



Airfoil fooMtry 



Figure 4.1.9 Airfoil Geometry and Measured and Predicted 
Loadings on the Pratt &. Whitney Tbansonic Cas ca d e. 


4.1.6 Energy Efficient Engine (E 3 ) Tirbine Lightweight Cascade 

The ability of a 2D version of the Rai code [TOMCAn - Domey, et al. (1992)] to compute flow 
for an airfoil over a range of incidence angles was verified by comparing its predictions against dam 
obtained by Sharma et al. (1982). Results from this verification effort, conducted under NAVAIR 
Contract #N00140-88-C-0677 at UTRC, are shown here to demonstrate the predictive capabilities 
of the code. Grids used in these calculations, containing 101*71 points for the inner grid and 41*21 
points for the outer grid are shown in Figure 4.1.10. Measured airfoil surface static pressure 
distributions compared to theoretical predictions from TOMCAT-2 and from VISCAS steady 
Navier-Stokes developed by Davis et al. (1986) in Figure 4.1.11 (a through f) over a range of 
jtyMfrnrp. angles. Results indicate excellent agreement between data and predictions. Predicted total 
pressure losses for the airfoil are compared to the experimental data in Figure 4 . 1 . 12 . Overestimates 
of losses at the design and the negative incidence angles are due to the assumption of fully turbulent 
flow in the analysis whereas in the experiment airfoil had large regions of laminar and traditional 
flows; turbulent flow at positive incidence angles is appropriate since leading edge separation bubbles 
induce transition in separation regions. 
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Comparisons of experimental data with theoretical predictions shown in this subsection clearly 
demonstrate that this code provides accurate estimates of airfoil loadings, losses and heat loads. The 
2D predictive capabilities of the code are, therefore, verified. 



a: Outer H—grid a: Inner O—grid 

Figure 4.1.10 Viscous Computational Grid for E 3 Lightweight TUrbine Blade 



Figure 4.1.11a Pressure Distribution for E 3 TUrbine Blade, 10* Incidence 



Figure 4.1.11b Pressure Distribution for E 3 liirbine Blade, 5 * Incidence 



Figure 4.1.11c Pressure Distribution forE 3 liirbine Blade, O' Incidence 
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Figure 4.1. lie Pressure Distribution for E 3 TUrbine Blade, —10’ Incidence 
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Figure 4.1.11 f Pressure Distribution for E 3 TUrbine Blade, —15 Incidence 



Figure 4.1.12 Total Pressure Loss Map for E 3 Lightweight Turbine Blade 

42. Verification of the 3D Steady Aspect of the Code 

Predictive capabilities of the code in a three-dimensional steady flow environment were 
verified by comparing theoretical predictions with benchmark quality experimental dataobtamed for 
an annular and a linear cascade. The emphasis of this effort '^^demonstrate that the ^ e P r f^ es 
reliable estimates of loadings, secondary flows and heat loads both on airfoil surfaces and endwalls. 
Studies were conducted to quantify the number of grid points on the pertinent features of 
three-dimensional flowfields in cascades. Limitations of computer resources on the NASA MSFC 
CRAY— XMP computer did not permit establishment of a gnd independent solution, furth 
evaluation of the impact on the grid numbers on the flowfield was, therefore, conducted as a part of 
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the NAVAIR Contract #N00140-88-C-0677 at UTRC using a version of the Rai code (Dorney, 
et al. (1992)) and one of the linear cascade geometries used in the present program. Highlights of the 
results from the NAVAIR contract are also discussed to demonstrate the predictive capabilities of the 
code. The imp*** of the number of iterations on the predicted flowfield was also investigated to 
identify degree of convergence needed to establish ultimate levels of losses and heat loads. A detailed 
discussion on the 3D code verification effort conducted under the present program was presented by 
Griffin and Belford (1990) and Griffin, Belford, Sharma and Ni (1991) in international conferences. 

42.1 Dring’s Annular Cascade - UTRC LSRR First Stator 

Detailed aerodynamic data consisting of spanwise airfoil loading distributions, midspan heat 
transfer coefficients, and flowfield information downstream of the stator are available for the UTRC 
LSRR first stator. Comparisons between predictions and experimental data are made to verify the 
heat tr ans fer, secondary flows, and performance prediction capabilities of the code. The code was run 
in a 3-D annular mode, and calculations were performed assuming a transitional flow over the airfoil 
surfaces. 

Airfoil Pressure Distributions . Comparisons of the calculated and experimentally measured 
airfoil pressure coefficients are shown in Figure 4.2.1 for the 2%, 50%, and 98% spanwise locations. 
The pressure coefficient is defined as 


C - P ~ P » 

p 1/2 qo0> j 

The agreement between the measured and computed results is excellent and is consistent with the 
results reported by Rai (1987). 
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Figure 4.2.1 Airfoil Pressure Distribution for UTRC LSRR First Stator for the 

(a) 2%, (b) 50%, and (c) 98% Spanwise Locations. 

Secondary Flows 

Secondary flow structure can be seen in plots of exit total pressure contours. The vortex 
structures contain fluid from both the endwall and airfoil boundary layers, and consequently represent 
regions of low total pressure. Exit total pressure contours for the LSRR first- stage stator are shown 
in Figure 4.2.2. The total pressure was measured and computed at a location of 17% of the airfoil 
chord aft of the stator trailing edge. Figure 4.2.2a shows the experimentally measur^ntours. F^ure 
4 2.2b shows total pressure contours calculated using a fairly coarse outer grid (the outer gr 
dimensions were 50 x 31 x 25) and an inlet total pressure that was constant across the span. The 
predictedresultsshowqualitative features of the measured data. Both show thepassage vortices, and 
die migration of these vortices toward the midspan. However, the computed wrticesj ire closer to 
their respective endwalls than those that were measured. The measured and predicted low total 
pressure regions compare well in terms of local loss. The maximum local loss of the dp secondly flow 
was measured to be Cm = -1.5 and predicted to be -1.7. The maximum local loss of the hub 
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secondary flow was measured at — 2-5 and computed to be — 2.1. The calculation yielding the contours 
in Figure 4 , ? ? ^ used the same grid as in the previous calculation, but a measured (spanwise - varying) 
inlet total pressure was imposed. The inlet boundary layer was measured to be 83% of the span at 
the hub and 11.67% at the casing. The locations of the computed and experimentally measured hub 
passage vortices are now in closer agreement The calculated casing passage vortex, although it has 
migrated further midspan, is still nearer the casing endwall than was measured. The shapes of the 
ryjniintf/i passag e vortices now more closely resemble the shapes of the measured vortices. The 
maximum local loss of the hub secondary flow was — 2.2. Additional spanwise planes were added to 
the grid in an effort to better resolve the endwall effects (as was done by Madavan ct al., 1991). Figure 
4.23d shows the exit total pressure contours calculated using the refined grid (50 x 31 x 31) and the 
measured inlet total pressure profile. This calculation produced results with much greater resolution. 
The maximum local loss of the tip secondary flow was still predicted to be —1.7, and the predicted 
maxim um local loss of the hub secondary flow was -23. The casing endwall vortex computed with the 
relatively fine grid is located further from the casing than its coarse grid counterpart, but nearer the 
racing than the measured vortex. Additional refinement in the spanwise direction is expected to 
further improve the calculated location of the casing vortex. However, due to computer memory 
limitations, additional spanwise refinement was not attempted during this study. One possible 
explanation for the discrepancy between prediction and measurement at the casing is the effect of 
concave curvature. Extra rates of strain introduced by normal pressure gradients, such as those due 
to surface curvature, have been shown to have large effects on viscous boundary layer development 
(Sharma and Graziani, 1982). The effect of surface curvature was not modeled in these calculations. 


Cuing 



Figure 4.22 Total Pressure Contours Downstream of the UTRC LSRR Fust Stator. 

(a) Experimentally Measured (Sharma, et aL, 1985), 

(b) Coarse Grid/Constant Inlet Total Pressure Calculation, 

(c) Coarse Grid/Measured Inlet Total Pressure Calculation, and 

(d) Spanwise- Refined Grid/Measured Inlet Total Pressure Calculation 
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Total Pressure Loss 

A typical measure of performance for internal flows is the total pressure loss through the 
passage. Comparisons of computed and measured total pressure loss through the LSRR first— stage 
stator are presented in Figure 433. Tbtal pressure loss coefficient is defined as 

r Ftp ~ Pi» 

l/2coO> s * 

where is drcumferentially averaged at each spanwise plane. The stator “exit” total pressure was 

computed and measured at a location 17% aft of the stator trailing edge. Figure 4.23a shows total 
pressure loss calculated using the relatively coarse grid that was di scussed in the secondary flows 
section. A total pressure that was constant across the span was imposed at the inlet Corresponding 
experimental data is also shown. Features represented by the measured loss data are predicted by 
RAI3DC. Both the experimental and computed results display the maxima and minima associated 
with the endwall secondary flows, but the computed and measured magnitudes and locations of these 
features do not agree very well. That this would be the case was highly predictable after observing the 
trends of the computed exit total pressure contours. For the coarse grid/constant inlet total pressure 
calculation, both the hub and casing passage vortices are positioned closer to their respective endwalls 
than the experimentally indicated location. The measured casing boundary layer is much thicker than 
the computed boundary layer. The calculation and experimental data agree reasonably well in the 
midspan region (the pressure loss coefficient at midspan was measured to be 0.14 and predicted to 
be 0.17). In an attempt to more accurately model the stator boundary layers and reduce the artificial 
loss generated in the grid overlap region, the outer grid was refined in the circumferential direction 
as was suggested by Rangwalla (1989). For a two-dimensional calculation, Rangwalla increased the 
number of circumferential grid lines from 31 to 71 and produced a midspan total pressure loss 
coefficient of 0.13. Unfortunately, due to limitations in computer memory in this study, it was not 
feasible to increase the number of circumferential planes by more than 120% as was done by 
Rangwalla. Grids containing a 20%, 40%, and 60% increase in the number of circumferential planes 
were generated. However, calculations using each of these grids yielded negligible differences in 
pressure loss. Apparently, the grid must be greatly refined in the circumferential direction in order 
to achieve any benefit in total pressure loss prediction. 

Predicted total pressure loss in Figure 433b was calculated using the coarse grid and the 
measured inlet total pressure profile. The magnitude of the computed and measured maxima and 
minima are now in dose agreement (except at 20% span where the code fails to predict the low total 
pressure loss for all cases and actually predicts a decrease in total pressure loss from 20% span to 
midspan). The predicted location of the relatively high loss region assodated with the hub passage 
vortex more dosely matches the data than the results of constant inlet total pressure case. The 
predicted location of the loss region coindding with the casing vortex moved slightly toward the 
mititpan The computed casing boundary layer appears to be slightly thinner than in the previous case. 
Figure 433c shows pressure loss coefficients calculated wit the spanwise— refined grid and measured 
inlet total pressure profile. The computed and measured magnitude and location of the relatively high 
loss region assodated with the hub vortex are not in excellent agreement. The predicted location of 
the loss region corresponding to the casing endwall vortex, though it has moved further from the 
raring , still exhi bited poor agreement with the measured location. The reason for this poor agreement 
is thought to be due to the effects of surface curvature, as was discussed in the secondary flows section. 

It must be mentioned that predicting total pressure loss required a relatively large amount of 
computer time. Although the airfoil loading predictions were converged after the residuals had 
dropped three orders of magnitude, and the Stanton numbers were converged after a residual drop 
of four orders of magnitude, the total pressure loss calculation required a drop in residuals of six 
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orders of magnitude. In addition, when a spanwise— varying inlet total pressure profile was imposed, 
a very small tim e step (one— fifth the size of the time step used in a constant inlet total pressure 
calculation) had to be employed to enable the solution to converge. Naturally, using the small time 
step increased the run time. 




(a) (b) 



Figure 4.23 Total Pressure Loss Through the UTRC LSRR First Stator. 

(a) Coarse Grid/Constant Inlet Total Pressure Calculation, 

(b) Coarse Grid/Measured Inlet Total Pressure Calculation, and 

(c) Spanwise -Refined Grid/Measured Inlet Total Pressure Calculation. 

Thermal Results 

A comparison between measured and computed midspan Stanton numbers is shown in Figure 
A2A. The Stanton number is defined as 


St - 


q" 
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Agr eement between the prediction and experimental data is excellent. 
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Figure 4.2.4 Stanton Numbers for the Midspan of the UTRC LSRR First Stator. 

Stanton number contours for the airfoil surface are presented in Figure 4.2.5. The pressure 
surface shows contour lines typical of a 2— D flow. The contours near the endwalls are attributed to 
the assumption of fully turbulent flow within 20% of the endwalls and not to secondary flows. The 
suction surface heat transfer, on the other hand, exhibits effects of a 3-D flow except for a region on 
either side of midspan which is primarily two-dimensional. Relatively high values of Stanton number 
and a steep gradient are seen at the leading edge and transition location on the suction surface. 
Additional regions of high heat transfer are located downstream of the leading edge near the endwalls. 
Avery steepgradient is noticed at that location. This region coincides with the 3-D flow of the passage 
vortices as they move from the endwalls onto the airfoil suction surface. 



Figure 4.15 Predicted Stanton Number (St x 10 3 ) Contours for the UTRC 
LSRR First Stator Surface. 

(a) Pressure Surface, 

(b) Suction Surface. 
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Stanton numbers for the endwalls are shown in Figure 4.2.6. The Stanton number contours 
upstream of the stator are basically parallel to the airfoil leading edge. This pattern implies the 
approach of a 2-D boundary layer. An area of high heat transfer is located at the stator leading edge 
where the passage vortices cause a rapid transport of fluid perpendicular to the endwalls. 



Figure 4.2.6 Predicted Stanton Number (St x 10^) Contours for the UTRC 
LSRR First Stator Endwalls. 

(a) Hub, 

(b) Casing. 


4J22 Langston Cascade 

Langston cascade represents one of the best documented benchmark quality data for an airfoil 
row representative of a turbine environment. Detailed aerodynamic experimental data were acquired 
for this cascade by Langston et al. (1977) at the UTRC in a large scale low speed cascade tunnel with 
three incoming profiles consisting of thin (1.8% of the span), nominal (13.6% of the span) and thick 
(50% of the span) boundary layers. Surface heat transfer coefficient data were subsequently acquired 
on this airfoil by Graziani et al. (1980) for the thin and the nominal inlet boundary layers. This cascade 
was tested with an aspect ratio of 1.08 and a pitch to chord ratio of 0.955. Test conditions involved 
ambient air entering the cascade at approximately 98 ft/sec. The Reynolds number based on axial 
chord (11.08 indies) and inlet velocity was 565,000. Additional details are given by Langston and 
Graziani. This test case was chosen to verify the loading, secondary flow and heat load prediction 
capabilities of the code. A grid refinement study was conducted to determine the sensitivity of 
predictions to grid points in the spanwise direction. The code was run with grids having ranges of 31 
to 49 points in the spanwise direction for the thin inlet boundary layer. A line of symmetry condition 
was imposed in the code for the nominal inlet boundary layer case and the code was run with girds 
having ranges of 16 to 39 points up to the midspan. The calculations were conducted both in a 
transitional and a fully turbulent mode. Most of the results discussed below were obtained with the 
finest grid and with calculations conducted in a transitional mode. 
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Surface Static Pressures and Surface Streamlines 

A comparison of measured and predicted airfoil surface static pressures is shown in Figures 4.2.7 
and 4.2.8 for the thin and the nominal inlet boundary layer. Excellent agreement is shown between 
experimental data and predictions. Figures 4.2.9 and 4.2.10 show a comparison of the measured and 
predicted surface static pressure coefficients at the endwall for the thin and the nominal inlet 
boundary layers respectively. Excellent agreement is shown between the data and predictions for the 
thin inlet boundary layer case. For the nominal inlet boundary layer, however, predicted results show 
larger impact of secondary flows than experimental measurements indicating an overestimation of 
secondary flows. 



Figure 4.27 RAI3DC Code Yields Excellent Agreement with the Airfoil 
Loadings at Various Spanwise Locations Measured by 
Graziani et aL (1980) for the Thin Incoming Boundary Layer 



I/ta 

Figure 4.28 RAI3DC Code Yields Good Agreement with the Airfoil 
Loadings at Various Spanwise Locations Measured by Langston 

etaL (1977) for the Nominal Incoming Boundary Layer 
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Ixparloantal Data UI3DC Calculation 



Figure 4.19 EndwaUCp Contours Thin Boundary Layer 


Experimental Data IH3DC Calculation 
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Figure 4.110 Endwall Cp Contours Nominal Boundary Layer 
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Figures 4.2.11 and 4.2.12 show a comparison of the measured and predicted streaklines on the 
airfoil suction sides for the thin and the nominal boundary layer cases. This comparison also indicates 
that whereas the thin boundary layer case is fairly well predicted; secondary flows for the nominal 
boundary layer case are overestimated which results in larger spanwise penetration of the separation 
line on the airfoil suction surface. 
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Figure 4.211 Theoretical Predictions Overestimate the Penetration of Separation 

Line on the Airfoil Suction Surface. A Possible Solution May Be to 
Increase the Number of Grid Points in the Spanwise Direction. 
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Figure 4.212 Penetration of the Separation Line on the Airfoil Suction Surface 
Overpredicted by the RAI3DC Cascade Code. Increasing the Number 
of Grid Points in the Spanwise Direction has a Favorable Effect on the 
Penetration Height at the Trailing Edge of the Airfoil Suction Side. 

Secondary Flows and Total Pressure Loss 

A comparison of the measured and computed total pressure loss coefficient at 10% axial chord 
downstream of the trailing edge of the cascade is shown in Figures 42.13 and 4.2.14 for the thin and 
the nominal inlet boundary layers. Excellent agreement is demonstrated for the thin inlet boundary 
layer case for the shapes and levels of pressure contours. Predicted loss contours for the nonunal 
boundary layer case, although they show similar behavior as the experimental data, indicate a large 
wake at the midspan and higher maximum total pressure loss coefficient than the experimental data. 
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MIDSPAN 



Figure 4.2.13 Total Pressure Loss Contours Downstream of Langston ’s Cascade 

for Thin Inlet Boundary Layer. 

(a) Experimental (Grariani et aL, 1980), 

(b) Predicted. 

MIDSPAN 



(a) EMDUALL (b) 


Figure 4.214 Total Pressure Loss Contours Downstream of Langston’s Cascade, 

Nominal Inlet Boundary Layer. 

(a) Experimental (Langston et aL (1977)), 

(b) Prediction. 

Predicted gap-averaged total pressure loss coefficients for the two inlet boundary layer cases 
are compared to the experimental data in Figures 4.2.15 and 4.2.16. Overall agreement between data 
and predictions is good except at the midspan for the nominal inlet boundary layer case. The 
agreement between data and predictions should improve if additional grid points are used in the inner 
grid and the effect of dilation induced by the secondary flow is accounted for through improvements 
in the turbulence model as suggested by Sharma and Graziani (1982). Limitations of computer 
storage available at the NASA MSFC CRAY-XMP computer did not permit further refinement of 
the grid especially in the airfoil surface normal direction. The effect of grid refinement on the 
predicted results was, however, conducted as a part of the NAVAIR contract #N00140— 88— C— 0677 
at UTRC by Doraey et al. (1982) by using a version of the Rai code for the thin inlet boundary layer 
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case by utilizing the NAS super computers and large storage super mini workstations. Results from 
Donley’s work, shown in Figure 4.2.17, indicate that almost 30% further increase in grid points are 
needed to provide an accurate estimate of profile and secondary losses in cascades. Significant 
im pr ovements in computer resources, especially by using super workstations, have, been achieved 
ri nw the present computations were conducted. Execution of reliable loss computations for turbine 
cascades is now within the reach of design engineers. It should be pointed out here that the code gives 
reliable results, provided the grid resolution is adequate. 



Figure 4.215 Total Pressure Loss, Gap Averaged, Langston’s Cascade 



Figure 4.216 Total Pressure Loss, GapAveraged, Through Langston’s 
Cascade, Nominal Inlet Boundary Layer 
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Figure4.217 Grid Refinement Study Conducted by DomeyetaL, (1992) 

Indicates Almost 550,000 Grid Points Are Needed to 
Accurately Resolve Losses in TUrbine Cascades 

Surface Heat Transfer Coefficients 

Measured and predicted heat transfer coefficients at the midspan of the airfoil for the thin and 
the nominal boundary layer cases are compared in Figure 4.2.18. The agreement between 
experimental data and predictions is fairly good; this agreement should improve with an increased 
number of grid points in the surface normal direction as pointed out above and as demonstrated by 
Domey et al. (1992). 

Distributions of Stanton number measured over the airfoil surfaces and endwalls for the two 
boundary layer cases are shown in Figures 4.2.19 through 4.2.22. Qualitative agreement is shown 
between the experimental data and theoretical predictions. Improvements in turbulence/transition 
models and increased number of grid points should further improve the agreement between data and 
predictions. 

Comparison of experimental data with theoretical predictions shown in this section clearly 
demonstrates the predictive capabilities of the Rai code. Lack of good agreement between the 
experimental data and theoretical prediction can invariably be attributed to the limitations of grid 
points used in the present effort. For engineering design execution, however, this code, even with 
limited grid resolution, provides excellent results which should allow use of the code in the design 
execution process for optimizing airfoil rows in rocket propulsion systems. 
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Figure 4.118 Midspan Stanton Number 
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Figure 4.2.19 Blade Surface Stanton Number Contours 
Thin Boundary Layer 
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Figure 4.220 Blade Surface Stanton Number Contours 
Nominal Boundary Layer 
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Figure 4.122 EndwaU Stanton Number Contours 
Nominal Boundary Layer 
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43 Tip Leakage Prediction Aspect of the Code 

Although initial version of the Rai code (ROTOR3) allows computations in the tip clearance 
region, these calculations did not provide stable results as noted by Rai (1987). Tip clearance levels 
used in computations by Madavan et al. (1989) were lower and by Dorney et al. (1992) were larger 
than the experimental results to yield stable overall solution. A study was conducted in the present 
program to verify the tip leakage flow prediction capabilities of the code. 

Flow was simulated in the Langston cascade with tip clearance flows on one of the endwalls. 
Ex p erimen tal dat a was available from this configuration from detailed measurements conducted by 
Dishart and Moore (1989). Computations were conducted with 39 points in the spanwise direction; 
which included 9 grid points in the tip clearance region containing 2.1% of the span. Initial calculations 
indicated insufficient grid resolution to resolve the overall flowfield but limitation of available 
computer storage space, on the NASA MSFC CRAY - XMP computer, did not permit inclusion of 
additional grid points. The calculations were, therefore, run to convergence to establish the level of 
accuracy achieved for the tip flows. 

A comparison of measured and predicted streaklines near the tip endwall for this configuration 
is shown in Figure 4.3.1. The predicted results yield fairly accurate description of the flow behavior 
near the endwalls. Separation and reattachment streaklines are fairly accurately predicted by the 
code; migration of flow from the pressure to the suction side is also well predicted. Measured total 
pressure loss contours 40% of the axial chord downstream of the cascade are compared to the present 
theoretical predictions in Figure 4 32 . The predicted results indicate a stronger endwall passage 
vortices than measured results, these results are similar to the ones discussed in the previous 
subsection for the nominal inlet boundary layer. The predicted secondary flow structures tend to 
merge into a high loss region which affects the convection of the tip leakage vortex. Distinct tip leakage 
vortex is, however, predicted as indicated in Figure 4.3.2 at the 10% and the 40% axial chord locations. 
The predicted tip clearance vortex is smaller than the experimental data primarily due to the influence 
of larger t han measured secondary flow structures. This exercise indicates that the code with a larger 
number of grid points in the spanwise direction will provide accurate simulation of tip clearance flows. 

It should be pointed out here that the tip clearance flows are dominated by inviscid migration 
of flows where the viscosity plays only a limited role. This point is illustrated in Figure 433 where the 
flow in a low turning cascade (Yaras et al. (1989)) was simulated by Staubach (1990) by utilizing an 
Euler code with 51 grid points in the spanwise direction. Good agreement between the measured data 
and Euler predictions demonstrates the need to use an increased number of grid points for tip 
clearance simulation; the need to account for viscous effects is not as important as originally 
anticipated. 
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Figure 4.3.3 3D Euler Calculation with 51 Grid Points in Spanwise 
Direction Can Resolve Tip Leakage Vortex ( Staubach (1990)). 

44 Unsteady Flow Prediction Aspect of the Code 

In addition to predicting the time- averaged loadings, surface streaklines and total pressure 
losses, an unsteady code must provide predictions of the following four well established unsteady flow 
features to ensure its proper verification: 

• unsteady pressure amplitudes on airfoils 

• segregation of hot and cold air in rotor passages 

• periodic elimination of rotor secondary flow vortices 

• unsteady variation in the transition point on downstream row airfoils and 
unsteady loss variation 

A number of simulations have been conducted in literature by using the code used in the present 
program to demonstrate its predictive capabilities. Results from some of these simulations are 
discussed below to provide evidence of the code verification. 

Rai (1987, 1989), Rai and Madavan (1988) and Madavan et al (1991) have provided sufficient 
comparisons with experimental data to demonstrate that the code provides accurate estimates of 
tim e— averaged flows for a turbine stage. Most detailed results for the UIRC LSRR were obtained 
by Madavan et al (1991) by using almost 1.43 million grid points to simulate three-dimensional flow 
through 3 vane and 4 blade passages; the experimental rig contains 22 vanes and 28 blades. Madavan 
et al. (1991) compared their results with those obtained by Madavan et al. (1989), which used 0.41 
million grid points and utilized 1 vane and 1 blade to simulate the flow through the same rig. Results 
from these publications, shown in Figures 4.4.1 through 4.4.4, clearly demonstrate that the code, with 
multi-passage simulation, provides very accurate predictions for: 

• spanwise distribution of time— averaged airfoil loadings (Figure 4.4.1), 

• amplitudes of unsteady pressures at midspan (Figure 4.4.2), 
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• time-averaged streaklines on airfoil surfaces (Figure 4.43), and 

• tim e— averaged total pressure loss contours at the exit of the vane and the blade 
(Figure 4.4.4). 




Figure 4.4.1(a) Spanwise Variation of Time-Averaged Pressure 
Distributions on the Stator. 

(a) 2% Span; 

(b) 25% Span; 

(c) 75% Span; 

(d) 98% Span. 

(Madavan et aL (1989, 91) Computations) 
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Figure 4.4.1(b) Spanwise Variation of Time-Averaged Pressure 
Distributions on the Rotor. 

(a) 2% Span; 

(b) 25% Span; 

(c) 75% Span; 

(d) 98% Span. 
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Figure 4.4.2 Pressure Amplitude Distribution on the Rotor at 
Midspan ( Madavan et aL (1989, 91) Computations) 
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Figure 4.4.3(a) Stator Surface Flow Visualization. Time-Averaged 
Limiting Streamlines From 

(a) Single-Passage Computations; and 

(b) Multi-Passage Computations. 

(Madavan et aL (1991)) 
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Figure 4.4.3(b) Rotor Surface Flow Visualization. Time-Averaged 
Limiting Streamlines From 

(a) Single -Passage Computations; 

(b) Multi-Passage Computations; and 

(c) Experimental Results. 

(Madavan et aL (1991)) 
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Figure 4.4.4 Relative Total Pressure Contours at the Exit to the 
Stator and the Rotor. 

(a) Single -Passage Computations; 

(b) Multi-Passage Computations; and 

(c) Experimental Results. 

(Madavan et aL (1989, 91) Computations) 


The multipassage simulation, conducted with 51 spanwisc grid points and with correct tip 
clearance levels also captures the tip leakage vortex structure which was not previously captured with 
a smaller tip clearance level by Rai. 


Pressure sides of turbine rotors operating in the presence of circumferentially inlet temperature 
profiles tend to have higher surface temperatures on the airfoil pressure sides than the 
time-averaged temperature at the rotor inlet. This flow phenomena, experimentally documented by 
Butler et al (1986), Sharma et al (1990) and Roback and Dring (1991), illustrates that the hot and cold 
stream of fluid tend to segregate in a turbine rotor. A number of unsteady numerical simulations (Rai 
and Dring (1990), Ni et al (1988), Ni and Sharma (1988), Krouthen and Giles (1988), Dorney etal 
(1991 1992), Thkahashi and Ni (1990, 1991) have been conducted to predict this flow phenomena. One 
of the most accurate simulations of this effect was conducted by Dorney, Davis and Sharma (1991) 
by using a 2D version of the code used in the present program. Results from this publication clearly 
demonstrate that segregation of hot and cold fluid in turbine rotors is an unsteady two-dimensional 
phenomena and it is well predicted by the code used in the present program as demonstrated by 
excellent agreement between experimental data and theoretical predictions in Figure 4.4.5. 


Of POOR Q.J.. , TY 
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Figure 4.4.5 A Radially Uniform Incoming TWo— Dimensional Hot Streak 
Yields Higher Surface Temperature on the Pressure Side and Lower 
Temperatures on the Suction Side (DomeyetaL (1991)) Indicating 
Segregation of the Hot and Cold Air in liubine Rotors. 

Periodic e limina tion of the rotor secondary flow vortices due to the interaction between the 
upstream vane and the downstream blade flowfield, experimentally documented by Sharma et al. 
(1983), was initially envisioned to be viscous effect. Unsteady, three-dimensional flow simulations 
conducted by an inviscid Euler code, however, were found to be sufficient to predict this flow feature 
as indicated in Figure 4.4.6 from Ni and Sharma (1989). It can, therefore, be safely concluded that an 
unsteady Navier Stokes code (like the one used in the present program) will predict this phenomena. 
Further expensive simulations are not needed to demonstrate this prediction capability. Review of 
existing numerical results from Madavan et al (1991) should dearly illustrate this flow behavior. 
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Figure 4.4.6 
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Relative Total Pressure at Exit of the Rotor from Both the Unsteady 
Simulation and Experimental Data Indicate That the Unsteady 
Euler Code (Ni & Sharma (1990)) is Sufficient to Predict 
Experimental Behavior. 


Periodic variation of the onset and the extent of transition has been observed by Pfeil & Herbst 
(1978), Doorley and Oldfield (1984) for an unsteady flow simulation in a linear cascade and by Hodson 
(1983) in a turbine rotor. This phenomena is most pronounced in moderate and low Reynolds number 
flow situations in the turbine and it should have little effect on the rocket turbines where Reynolds 
numbers are very high. An unsteady transition model has been implemented in the pr esent code as 
dipleg**! in Section 3 of the present report. Flow simulations were conducted for the UTRC LSRR 
by using this transition model. Large leading edge overspced for the rotor, however, always tripped 
the boundary layer to a turbulent state. This effect, therefore, cannot be illustrated in this subsection 
of the report. Code application results shown in Section 5.2 do, however, discuss the impact of this 
phenomena on the overall unsteady and time-averaged losses. 

TVvo— and three-dimensional flow simulations were conducted in the present program for the 
UTRC LSRR rotor. Measured flow properties from the upstream stator (vane) were specified as inlet 
boundary conditions. Results of the two-dimensional calculations are similar to those obtained by 
Rai (1987) and Madavan et al. (1991) and some of these results are shown in Figures 4.4.7 and 4.4.8. 
Three-dimensional flow calculations were also conducted for the rotor by using a uniform upstream 
and measured upstream boundary conditions from the stator exit. Limitations of the available 
computer storage on the NASA MSFC XMP contained the spanwise grid density to 39. In addition, 
a fully converged solution has not yet been achieved so the results are not discussed. Even converged 
results from these simulations are not expected to be any superior to those obtained by Madavan et 
al (1991) because of fewer number of grid points used in the present simulation. It should, however, 
be pointed out that verification of the unsteady flow predictive capabilities of the code has been 
demonstrated by results shown in this section. 

In summary, results shown in this section clearly demonstrate the predictive capabilities of the 
code used in this program. The code provides accurate estimates of losses, heat loads, loadings, 
unsteady pressure amplitudes, secondary flows, tip leakage vortices and segregation of the hot and 
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cold air in turbine passages. The only requirement for accurate predictions is proper grid resolution 
and computer resources to ensure convergence of calculations. 


Static Pressure 


Pressure Coefficient 




tSfra 


New Model 


New Model 


Figure 4. 4. 7 Unsteady Two-Dimensional Computations for the UTRC LSRR 

Rotor Indicating Unsteadiness in Static Pressure, Skin Friction and 
Losses Larger Effects Indicated Using Transitional Model 

Developed in the Present Program. 
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For stsady-stats solution 
(data is gap-avsragsd) 


For unsteady solution 
varies with pitch and radius 


Figure 4.4.8a Computational Grid and Inlet Profiles Used in Simulating Steady 

and Unsteady 3D Flows Through the UTRC Rotor. 





Figure 4.4.8b A Comparison of Measured and Predicted Time— Averaged Airfoil 

Surface Static Pressure Coefficients for the UTRC LSRR Rotor 
Indicating Good Agreement 
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5. CODE APPLICATION 


The verified CFD code was applied to simulate flow through two configurations to demonstrate 
that it can be used to compute the flow through both axial and radial turbomachines. This exercise 
indicates that the CFD code developed under the present program can be applied in the design 
optimization studies for rocket turbopumps to yield high performance and improved durability 
hardware. The results obtained from this exercise are discussed below. 

5.1 Railly’s Radial Impeller 

Rainy’s impeller (Ekerol & Railly (1982)) was selected as one of the configurations in the code 
application task of the program to show that the code will provide accurate estimates of flows in radial 
(centrifugal) machines. This impeller is a low speed radial cascade with an outer diameter of 45 indies 
as shown schematically in Figure 5.1.1. Air enters this impeller radially after passing through a 
vaneless inducer which directs the flow from the axial to the radial direction. Large contraction in the 
inducer ensures uniform flow at the impeller inlet. The pitch to chord ratio of the impeller is 0.78 and 
the Reynolds number, based on radial chord (7.72 inches) and an exit relative velocity of 433,000. The 
span of the impeller is constant (3.5 inches). The airfoils are shrouded to ensure absence of leakage 
flows. 


Rig schematics CompiMional grid 



Figure 5. 1. 1 Code Application — Railfy’s Impeller. 

Calculations were conducted in a fully turbulent mode by assuming a two-dimensional flow at 
the mid-span; effects of endwall boundary layers were accounted for through the use of stream-tube 
contraction ratio. Computational grids used in the code are shown in Figure 5.1.2. Predicted 
streaklines indicate separation of the flow on the airfoil suction side as indicated in Figure 5.13. 
Measured boundary layer thickness distribution along the airfoil suction side is shown in Figure 5.1.4. 
Experimental data indicate a large gradient in the boundary layer thickness around 66% of the chord 
on the airfoil suction side which is in close proximity of the onset of the separation zone predicted by 
the code (Figure 5.13) This comparison indicates that the code provides an accurate estimate for one 
of the most dominant features of the flowfield in radial impellers. Three-dimensional flow 
simulations were also initiated for this impeller; the solution was, however, not run to convergence 
due to the unavailability of computer resources. 
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Figure 5.1.2 


Grid Used in the Simulation. 
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Figure 5.1.3 Code Application - Railfy’s Impeller 
(2D Steady Calculation). 
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Figure 5.1.4 Code Application - Rally’s Impeller 
(2D Steady Calculation). 

52 Pratt & Whitney Full Scale Tbrbine 

This configuration was selected to demonstrate that the CFD code used in the present program 
provides a more accurate estimate of unsteady loads on airfoils than the Euler code; the latter code 
is extensively used in the current design process for turbomachines. 

The need to establish the relative accuracy of the unsteady loading prediction capability of an 
Euler and a Navier- Stokes code became apparent in the design of the Space Shuttle MainEngme 
(SSME) alternate turbopump development (ATD) high pressure oxygen turbopump (HPOT) and 
high pressure fuel turbopump (HPFTP) turbine blades. Analysis of the first stage blades of each 
turbine indicated possible resonance problems in crucial operating ranges of the turtopraps. 
Unsteady aerodynamic simulations were conducted for the two turbines (Griffin & Rowey(1993)) to 
support further investigations of the dynamic responses of the first stage turbine blades. These 
simulations were conducted at the mean radius of the two turbines by utdmng two dimensional 
unsteady Euler (Ni et al (1990)) and Navier Stokes (Gundy-Burlet et al(1991)) codes. Both axles 
were found to yield almost identical time-averaged loadings. Unsteady axial and tangential loadings 
for the first blades of the two turbines were, however, found to be very different as shown in Figure 
5.2.1, which also shows a comparison of the time-average loadings for the two first blades. The 
unsteady loads predicted by the Navier-Stokes code were almost an order of magnitude larger than 
those predicted by the Euler code. The main objective of the present application is to evaluate the 
unsteady loading prediction capabilities of the unsteady Euler and Navier Stokes codes for a subsonic 
axial flow turbine. 

A brief discussion of the experimental turbine rig, where appropriate data are available to verify 
the predictive capabilities of the two codes, is given below. 
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Figure 5.2.1 Numerical Simulations Conducted For The HPOT And HPFTP 
Thrbines (Griffin A Rowey (1993)) By Using Two-Dimensional 
Unsteady Euler And Navier— Stokes Codes Show: 

I) Almost Identical Time— averaged Loadings; But 

II) Smaller Unsteady Loads From The Euler Code Than Those 
From The Navier— Stokes Code. 


Experimental Rig 

The experimental rig consists of a two stage turbine operating at an overall pressure ratio of 4.5 
and a Reynolds number of about 700,000 for the first stage stator. The schematics of the rig and 
measurement locations are shown in Figure 5.2.2. Airfoil surface static pressures are measured at 
various spanwise locations for all airfoil rows. Overall two stage turbine performance is measured by 
traversing the flow at the exit of the turbine, using pressure and temperature rakes in four quadrants. 
i wading edge Kiel head probes at the second stator inlet are used to define the first stage performance. 
Unsteady total pressure data are acquired downstream of the turbine to quantify interaction between 
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airfoil rows. Laser Doppler Velocimetry (LDV) is used to acquire unsteady velocity data in between 
the airfoil rows. 



Figure 5.2.2 Schematics of the IWo— Stage TUrbine Rig Along With the 

Measurement Locations. 


Unsteady surface static pressures were measured at several locations at the mid-span of the 
second stator. The measurements were obtained using infinite tube probes coupled to conventional 
airfoil surface static pressure taps. The infinite tube probe technique permits unsteady measurements 
to be maHf on the surface of the airfoil without altering the surface contours, as would be required 
with the installation of surface mounted high-response pressure transducers. It also removes the 
transducers from direct contact with the flow field, thus providing a more durable high-response 
instrumentation system. 


In this application, the unsteady static pressures were obtained at six locations in a stator 
passage. Measurements were made at 15%, 50%, and 87% chord on the suction side of one airfoil and 
at 40%, 85% and 100% chord on the pressure side of an adjacent airfoil. The infinite tube probes were 
located outside the turbine case 16 to 20 inches from the pressure taps, resulting in a usable frequency 
response of up to 1 1 kHz at the rig operating pressures and temperatures. This is above the rotor blade 
passing frequencies of 6 to 7 kHz but not adequate to capture any higher harmonics. (Note: work is 
currently in progress to reduce the sense tube line lengths to improve the response of the probes to 
over 20 kHz). 


lb de termin e the pressure fluctuations at the airfoil surface measurement location, the data 
from the infinite tube probes must be compensated for the attenuation and phase shift m the 
connecting Unes. Typically, one hundred revolutions of data are acquired for each test point. Each 
revolution of data is then frequency compensated using the technique described by Nyland etal(1971 ). 
The compensated data is then ensemble averaged to produce the average periodic and random 
components of the data for one rotor revolution. Levels of periodic unsteadiness are obtained by 
examining the fluctuations of the periodic signals and levels of random unsteadiness are obtained by 
camming the mean levels of the random signals. 
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Solution Procedures (CFD Codes & Models) 

Euler Code 

Three dimensi onal multistage unsteady Euler flow solver, developed by Ni (1989,90), is used to 
co nduct in viscid numerical sim ulati ons. This solver utilizes an explicit time— accurate finite volume 
numerical scheme with the cell— vertex centered solution algorithm to integrate the unsteady Euler 
equ at i on s to provide a solution for flow in multistage turbomachines. The scheme is second order 
accurate both in time and space on a smooth computational mesh. Detailed numerical equations for 
im pl ementin g the scheme and methods for applying various types of boundary conditions are given 
inNi et al(1989). For time— accurate simulations of multistage flow, a cubic interpolation method is 
used for tr ansf er of flow information across the interface boundaries dividing stationary and rotating 
airfoil rows as discussed by Ni et al (1989,90). A surface drag force model similar to the one described 
by Denton(1990), is used to provide the effect due to viscosity near airfoil and endwall surfaces. The 
boundary conditions applied for the present simulations are flow tangency on solid surfaces, 
prescribed spanwise and tangential distributions of total pressure, total temperature and flow angles 
at inlet to the multistage turbine. The number of grid points used in simulations and their distributions 
along with the operating conditions are given in Figures 523 & 5.2.4 for the two— and 
three-dimensional simulations. 



Figure 5.23 Computational Grid Used to Simulate TWo— Dimensional Unsteady 

Flow Through the First Rotor and Second Stator of the Rig. 
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Figure 5.24 Computational Grid Used to Simulate Three-Dimensional Steady and 
Unsteady Flow Through the Multistage TWo— Stage Rig by Using the 

Euler Code. 

Navier— Stokes Code 

TWo— dimensional version of the three-dimensional unsteady Navier Stokes code, developed 
by Rai (1989), is used to conduct the viscous numerical simulation. The numerical procedure consists 
of a time marching, implicit, third order spatially accurate, upwind, finite difference scheme. The 
invisdd fluxes in the code are discretized according to the scheme developed by Osher (1982). The 
viscous fluxes are calculated using standard central differences. An alternative direction, 
approximate— factorization technique is used to compute the time rate change in the primary 
variables. In addition, an inner Newton iteration is used to increase stability and to reduce 
linearization errors. The code is modified to account for the ‘stream-tube’ (‘H-Ratio’) effects in a 
manner proposed Rangwalla et al (1991). 

The two- layer Baldwin-Lomax (1978) turbulence model is used to compute the turbulent 
viscosity. This turbulence model is modified to account for the transitional nature of the airfoil suction 
surface boundary layer by using the intermittency factor approach suggested by Sharma (1987). 
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Modified TUrbulent Viscosity = F * Tbrbulent Viscosity B— L 
F = Interim ttency Factor 

= l.-exp (-( R0**2J-R0c**2J))/R0c**2^8) 

R0 = Reynolds Number based on momentum loss thickness of the airfoil boundary layer 

R0c = Critical Reynolds Number - 40/Ili**0.5 

Hi = Free stream turbulence level 

Free stream turbulence level is calculated by using entropy value at the edge of the stator 
boundary layer. Circumferential variation of flow properties at the interface between the rotor and 
the second stator is used to determine a relation between entropy and turbulence for application in 
the second stator. TUrbulence level is specified as an input value to the rotor. 

Multiple zonal grids are used to discretize the rotor/stator flow field and to facilitate relative 
motion of the rotor (Rai(1989)). A combination of O— and H— grid sections are generated in the 
blade— to- blade direction extending upstream of the rotor leading edge to downstream of the stator 
trailing edge. Algebraically generated H- grids are used in the region upstream of the leading edge 
to downstream of the trailing edge and in the inter blade region. The O— grids, which are body fitted 
to the surfaces of the airfoils and generated by using elliptic solution procedure, are used to properly 
resolve the viscous flow in blade passages and to facilitate application of algebraic turbulence models. 
Computational lines within the O— grids are stretched in the blade— normal direction with a fine 
spacing at the wall Figure 5JZS illustrates the grid topology used in the present simulation. Boundary 
conditions used in the current simulation are no slip condition on solid surfaces; prescribed spanwise 
and tangential distributions of total pressure, radial and tangential flow angles; and Reimann 
invariant at the inlet; prescribed static pressure at the exit 



Figure 5.Z5 Computational Grid Used to Simulate Viscous Flow Through the 

Rotor and the Stator by Using an Unsteady Navier— Stokes Code. 
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Discussion of Results 

Unsteady two— dim ensional flow simulations were initially conducted to evaluate the steady and 
unsteady loading prediction capabilities of the Euler and the Nayier— Stokes codes. Predicted 
envelope of unsteady surface static pressures on the stator from these simulations are shown in Figure 
5.2.6. It is apparent from this figure that the Navier Stokes code predicts larger magnitudes of unsteady 
pressures than the Euler code over the entire airfoil surface. This result is similar to that obtained in 
AID turbines as shown in Figure 5.2.1. A comparison of the experimental data with predicted results 
in d j ratpg that even the Navier- Stokes code underestimates the unsteady pressures on the airfoil. In 
addition, analysis of numerical results indicates that the two-dimensional simulation of the turbine 
has yielded a lower inlet Mach number and a larger negative incidence angle for the stator than these 
were present in the experiment. 


o »teady data 
H UNSTEADY DATA (AMPLITUDE) 
SIMULATIONS 


20 EULER 20 NAVIER -STOKES 



Two-Dimensional Unsteady Flow Simulations Show Larger Levels 
of Unst ea din ess From The Navier— Stokes Code Than Those From 
The Euler Code. Both Simulations, However, Do Not Accurately 
Model The Inlet Condition To The Stator. 


A review of the steady three-dimensional multistage Euler flow simulations of the entire two 
stage turbine showed that the mid-span region of the second stator was strongly influenced by the 
endwall secondary flows and flowpath divergence. These steady simulations yielded excellent 
agreement with the time-averaged loadings on airfoil surfaces, predicted results at the mid-span of 
the first rotor and the second stator are compared to the experimental data in Figure 5.2.7. Unsteady 
three-dimensional flow simulations were then conducted for the two stage turbine using the Euler 
code. Predicted envelope of unsteady pressures for the second stator from this simulation is shown 
in Figure 5.2.8, amplitudes of measured unsteady pressures are also identified in this figure. Results 
shovel in this figure indicate that the 3D Euler code predicts larger amplitudes of imsteady pressures 
on the second stator than the 2D code primarily because it provides more realistic inlet flow conditions 
to the second stator. The predicted magnitudes of unsteady pressure are, however , still lower than the 
experimental data. 
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Figure 5.27 Three-Dimensional Steady Multistage Euler Code Yields Good 
Agreement with the Airfoil Loading Data for the Rotor and the Stator. 
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Figure 5.28 Three-Dimensional Multistage Unsteady Euler Code Predicts Higher 

Unsteadiness Than the 2D Code, But it Still Underestimates Unsteady 
Pressure Amplitudes on the Airfoil Suction Side. 

One is tempted to conclude at this stage that only a 3D unsteady Navi er Stokes code can provide 
predictions of unsteady pressures in this experiment. Computational resources required to 
e xecute a 3D unsteady Navier Stokes for this configuration are, however, prohibitive. The approach 
undertaken in the present investigation was to use the two-dimensional unsteady Navier Stokes code 
which was modified to account for the endwall flow effects through the use of 'steam tube’ (H— ratio) 
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variation. The distribution of stream tube variation was calculated from the numerical results 
obtained from the 3D steady multistage Euler flow simulation. 

Steady two-dimensional flow simulations were first obtained by using the cascade version of the 
unsteady Navier Stokes code for the rotor and the stator with the ‘H-ratio’ variation calculated above. 
Calculated airfoil surface static pressure distributions from these cascade simulations are coinPfnxi 
to the experimental data obtained by using steady pneumatic instrumentation m Figure 52.9. Oood 
agreement is shown between the experimental data and predictions. 



Figure 5.2.9 Steady Loadings on the Rotor and the Second Stator Well Predicted by 2D 

Steady Navier Stokes Code Modified to Account for “H— ratio” Effect. 


T\vo— dimensional unsteady flow simulations were then conducted for the rotor and the second 
stator by using the 2D Navier Stokes code modified for the ‘H-ratio’ effect. The predicted envelope 
of unsteady pressure from this simulation for the stator is compared to the experimental data in Figure 
5.2.10. Measured amplitudes of unsteady pressures are also identified in this figure. Good agreement 
between the unsteady experimental data and predictions for the stator clearly demonstrates that a 
Navier-Stokes code is needed to estimate unsteady loads in turbine rows. This result also shows that 
a two-dimensional code, modified to account for ‘H-ratio* effects, is a cost effective altema * lve *° 
three-dimensional unsteady Navier Stokes codes, at least for estimating the levels of unsteady loads 


on airfoil surfaces. 
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Figure 5.210 Envelope of Lootings on the Stator Predicted by Using 

2D Unsteady Navier— Stokes Code with “H— ratio” Modifications ; 
Amplitudes of Unsteady Pressures on the Airfoil Suction Side Fairly 

Well Predicted. 


Numerical results from the above unsteady flow simulation were processed to establish the levels 
of unsteady losses for the second stator as it was affected by the unsteady interaction with the upstream 
rotor. Periodic variation in the loss for the second stator is shown in Figure 52.11. The second stator 
is found to yield almost 50% higher time-averaged loss in the unsteady flow environment than the 
loss calculated for that airfoil in a steady flow environment. The predicted increase in the 
time— averaged loss for this airfoil is similar to that measured by Hodson (1983) as shown in Figure 
52.12. The predicted periodic variation in unsteady losses is similar to those experimentally measured 
for the second stator in the UTRC LSRR as shown in Figure 52.13. Results from this application 
indicate: 


• An unsteady Navier Stokes code provides a more accurate estimate of unsteady 
loads in a multistage turbine environment than an Unsteady Euler code. 

• 2D uns tead y Navier Stokes codes, modified to account for *H- ratio' effect, are 
a cost effective alternative to more expensive and computer intensive 3D 
unsteady Navier Stokes codes to estimate unsteady loads in turbines. 

• Time — averaged losses for airfoils in an unsteady flow environment are larger 
than estimated for those airfoils in a steady flow environment and unsteady 
Navier Stokes codes can model this flow phenomena accurately. 
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Figure 5.2*11 Two — Dint ensiofuil Unsteady Navier Stokes Code Predicts Periodic 

Variation in Loss for the Stator as it is Influenced by the Upstream Rotor. 
Tune-averaged Loss for the Stator is Almost 50% Larger Than 
Calculated for this Airfoil in a Steady Flow Environment 
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Figure 5 112 Measured Streamwise Distribution of Tune-Averaged Boundary Layer 
Momentum Loss Thickness (Hodson (1983)) Show Larger Values For 
Rotors (Unsteady Environment) Than Those Measured For The Same 
Airfoil Sections In a Steady Cascade Environment. The Rotor Data Are 
Bracketed By The Transitional and Fully TUrbulent Calculations (Sharma 
etaL (1988)). 
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Figure 5.2.13 Mid— Span Losses On The 2nd Stator As Influenced By The 
Unsteadiness Generated By The Upstream Rotor Airfoil Wake. 
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6. CONCLUSIONS AND FUTURE DIRECTION 


CFD code developed by Rai (1987) at NASA ARC has been modified to accountfor transmonal 
flows on airfoil surfaces. Improved turbulence models have been implemented in the ® od ® 
account of physical phenomena such as surface roughness, rotation, turbulence level and extra rates 
of strain on the development of viscous flows in turbopumps. This code has been shown to yield good 

estimates of: 


• Airfoil loadings, heat loads, losses and secondary flows in a two- JJ«*a 
three-dimensional steady flow environment (Griffin and Belford (1990), 
Griffin et al. (1990)). 


• Tip leakage flows. 

• Unsteady loads, time- averaged loadings and flowfield in an unsteady flow 
environment. 

• Unsteady variation of losses and increased time-averaged loss than measured 
for airfoils in a steady flow environment. 


• Flow separation in a radial impeller. 

Grid sensitivity studies conducted, in the present program and tho^OTnducted by Dorney et al. 
(1992) by using a version of the present CFD code, indicate that almost 550,000 grid points are n 
to accurately estimate profile and secondary losses in turbines. 

Studies conducted in the present program and those conducted by other investigators using 
versions of the present code can predict following known effects of unsteadiness m turbine stages: 


• Segregation of hot and cold air in turbine rotors (Domey, Davis and Sharma 
(1990)). 

• Periodic elimination of rotor secondary flow vortices in a turbine stage (Ni and 
Sharma (1990)). 

• Amplitudes of unsteady pressures in a turbine stage (Rai and Madavan ( 1989)). 
Present code yields more realistic estimates of unsteady loads than an unsteady 
Euler code (Sharma et al. (1992)). 


• Time-averaged losses are higher in an unsteady flow environment than 
measured for the same airfoil in a steady flow environment (Sharma et al. 

(1992)). 


Future work in this area needs to focus in numerical investigations to isolate and identifyloss 
generation mechanisms to provide guidance to designers. These investigation may focus on 

identifying: 


• Loss production mechanisms due to secondary flows in airfoil rows. 


• Loss production and control mechanisms due to tip leakage flows. 

• Loss production due to steady and unsteady flow interactions in multi-stage 
machines. 


Ill 
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Foreward and Summary 


This Document is a manual describing the use of a 3D Navier-Stokes flow 
solver, RAI3DS, which was developed under NASA Contract NAS8-36950, 

"Three Dimensional Turbopump Flowfield Analysis". The Program Managers 
of this contract are Lisa W. Griffin, of NASA/MSFC, and Dr. Om P. Sharma, of 
Pratt & Whitney Aircraft. This Document accompanies the final report of the 
contract, in which the methodology of the solution algorithm is presented. The 
focus of this manual is the mechanics of using the system, accompanied by 
some samples of jobstreams, input, and output. Some details of major modifications 
made to the original version of this program, ROTOR 3 (M. M. Rai, NASA/Ames, 
1987) are also included. 


1 


Overview of RAI3DS System 


The purpose of RAI3DS is to solve for the flow through a turbopump. It is the 
result of work in support of NASA Contract NAS8-36950, the main tasks of which 
were to verify, enhance, and apply the three-dimensional, unsteady Navier-Stokes 
flow solver ROTOR3, written by Man Mohan Rai, NAS A/Ames. ROTOR3 was 
designed for a specific experimental single-stage axial turbine, for which it gave good 
predictions of the experimental data (Ref. 1). In order to utilize the predicitve 
capability of ROTOR3 as a turbomachinery design tool, modifications were necessary 
to make the code applicable to general cases, and also to reduce the computational 
resources required for its execution. As a result of this work, RAI3DS will solve for 
one, two, or three airfoil rows, up to four airfoils in each row, in two or three 
dimensions, with or without tip clearance, planar or annular, steady-state or unsteady. 
It can also be applied to radial flow machines as well as axial. Additional features 
include capability for imposing a line of symmetry condition in three-dimensional 
axisymmetric cases; the option to model stream tube contraction effects in 
two-dimensional cases; and the capability for specifying an incoming wake and/or 
boundary layer. 

RAI3DS consists of three major modules: the grid generator, the flow 
initializer, and the flow solver, each having different forms of input and output. The 
grid generator and flow initializer run on the Iris workstation, and the solver runs on 
the Cray-XMP. In addition, there is a shell script, run on the Iris, that generates the 
UNICOS job stream needed to mn the flow solver on the Cray. Each of these 
modules are explained in detail, including sample input and output, in the following 
sections of this manual. 


1.1 



Grid Generator 


The solution scheme of the RAI3DS flow solver requires two overlaid grids for 
each airfoil. A dense inner grid surrounds the airfoil with lines normal to the surface, 
and is overlaid by a less-dense outer grid of horizontal and vertical lines extending to 
the boundaries of the solution domain. Between each airfoil, the outer grids are 
slightly overlaid at the interface. This grid scheme is explained in more detail by Rai 
(Ref. 2), but here is an example of a radial cross-section of a stage with rotation about 
the X-axis: 



In actuality, the outer H-grid extends throughout the solution domain; the innermost 
portion is not shown here so as not to obstruct the view of the inner O-grid. The grids 
are stored in Cartesian coordinates, with X-, Y-, and Z- coordinates for each node, 
represented by (IJ,K). For the outer grids, I varies with X and J varies with Y; for 
the inner grids, I varies with airfoil circumference and J varies with distance from the 
surface; K varies with radius,which is normal to the page. 

The procedure to generate these grids is as follows: for each airfoil, the outer 
H-grid is generated first, it’s boundaries in the Y-direction following the curve of the 
airfoil’s camber line, which is generated from the X-Y coordinates of the airfoil 
surface. The outer boundary of the inner O-grid is then generated, and the O-lines are 
distributed from the airfoil surface out to the boundary. The O-lines are connected 


with lines normal to the airfoil at the surface and curving out to the boundary. This is 
done at five spanwise locations, and then repeated according to the user-specified 
radial distribution of the grid. For planar cases, the grid is then unwrapped at each 
radial location so that each grid plane has a constant Z-value. For 2-dimensional 
cases, all grid planes except that at midspan are eliminated, and the midspan grid 
plane is copied twice, forming a grid of three identical planes at three different 
Z-locations. For radial cases, the planar grid (i.e., having constant Z-values at each 
X-Y plane) is transformed to curve around the Z-axis, as shown here: 



For three-dimensional cases with rotation, RAI3DS will solve for the flow over 
the tip of the rotor, requiring an extra grid for the clearance region. This grid is also 
an O-grid, having as its outer boundary the edge of the rotor tip and as its inner 
boundary an ”0" collapsed into a line. The O-lines are connected by lines that are 
normal to the airfoil surface in the X-Y plane, which are extensions of the normal 
inner grid lines: 



The above shows the first few lines of the inner O-grid meeting with the tip clearance 
grid at the airfoil surface. A tip clearance grid is generated at each grid plane above 


2.2 


the tip of the rotor, for which the radial distribution is specified by the user. 

For some unsteady cases, it is advantageous to use a small H-grid to represent 
the exit section of an airfoil row, rather than solve for the entire row. This is what is 
referred to as a "pseudo-airfoil" both in this document and in the RAI3DS system 
codes. The feature was added for the purpose of simulating first rotor interaction with 
the wake coming off of a first stator, without having to compute the flow through the 
first stator. A small, stationary H-grid was generated to represent the exit section of 
the first stator, and a wake was defined at the inlet and held as the boundary condition. 
The rotor grids were then moved past this small stationary grid just as they would if 
there had been a full stator grid upstream of them. This type of simulation saves on 
computation as well as storage requirements, but is only applicable to cases in which 
there are no unsteady effects present at the location of the inlet to the pseudo-airfoil 
grid. Here is an example of a pseudo-stator grid preceding a rotor grid: 



location of 
stator trailing edge 

The program requires the following input: the geometry of each airfoil, in X-Y 
coordinates, at several spanwise locations, the radius at which each of these is defined, 
the number of blades in each row of airfoils (i.e., the total number of airfoils in the 
disk), and the upstream and downstream boundaries of the solution domain of each 
airfoil, although these are subject to change during execution of the program. For 
each airfoil , this information is read in, and the user is prompted for information to 
tailor the grid to satisfy the requirements of the particular case. The user will have the 
opportunity to change grid density in all three directions, scale the airfoils relative to 
one another, and change upstream and downstream boundaries. Obtaining the 
optimum grid for a case generally requires several iterations of the program, even for 
those familiar with the effects of changing each variable. For this reason, the 
generated grid is output in a PLOT3D-readable file to allow the user to view the grid 


2.3 


and determine the necessary modifications, if any. 

To best explain the various options available, a sample run is presented with 
detailed explanations accompanying each user prompt. This sample is a 3D, planar 
rotor with tip clearance, its geometry being the same as the first rotor of the United 
Technologies Research Center Large Scale Rotating Rig (LSRR) (Ref. 3), and a 
pseudo-stator, meant to demonstrate nearly all possible options of the code. 

Obviously, much of this would be bypassed in a simpler case. The input file, the 
format of which is described at the end of this chapter, is expected by the program to 
be in FORTRAN unit 2. This is the only necessary file assignment, and the user may 
then invoke the program. 

% grd3ds 

ENTER 1 FOR A RADIAL IMPELLER CASE, RETURN FOR AXIAL 

HOW MANY FULL AIRFOIL ROWS ? 

1 

This is the number of full airfoil grids you want to generate, regardless of whether 
you want a pseudo-airfoil grid upstream. You must have coordinates for this 
number of airfoils in the input file. 

ENTER TYPE OF 1ST FULL AIRFOIL ROW, 1=STATOR, 2=ROTOR 

2 

"1st" means leftmost. 

PSEUDO-AIRFOIL UPSTREAM ? (0=N0, 1=YES) 

1 

2D OR 3D? (ENTER ’2’ OR ’3’) 

(NOTE 3D IMPLIES TIP CLEARANCE IF SOLVING FOR A ROTOR) 

3 

Note that if you want to solve for a stationary airfoil with tip clearance, refer to it 
as a rotor and define rotation speed to be zero (this is done later ). 

PLANAR CASE? (0=N0,1=YES) 

1 

This applies to 2 or 3 dimensions. 

ENTER NMBR OF AXIAL POINTS FOR PSEUDO-AIRFOIL GRID, 12 

10 

ENTER JMAX FOR PSEUDO-AIRFOIL GRID, 12 

28 

The spacing in the Y -direction in the pseudo grid should be as close as possible to 
that of the neighboring airfoil grid; this will depend on the pitch ratio of the rows. 
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ENTER 1 TO CHANGE AIRFOIL 2 K-PLANE MESH DIMENSIONS 
(INNER: 101,21 DSBOD = .0005 
OUTER: 50,31 
TIP: 101,11) 

1 

The above are defaults. 

ENTER IMAX (AXIAL PNTS) FOR OUTER GRID, FORMAT=l2 
(DEFAULT=50, MAX=61) 

ENTER 1 TO CHANGE IBEG AND IEND, OR RETURN 

1 

This allows changing the axial boundaries of the inner grid relative to the airfoil. 

ENTER % AX CHORD FOR IBEG AND IEND (CURRENTLY=.1 ) 

.08 

Change them from 10% bx to 8% bx, i.e. closer to the airfoil l.e. and t.e.. 

ENTER JMAX (CIRCUM. PNTS) FOR OUTER GRID, FORMAT=l2 
(DEFAULT=31, MAX=31) 

21 

This depends on JMAX for the pseudo grid, also pitch ratio; in this case, the 
pseudo grid will have 28 points, and there will be 3 pseudo-stator grids for 4 rotor 
grids, so 21 points are chosen so that the spacing will be the same for each row. 

ENTER 1 TO CHANGE JBEG AND JEND, OR RETURN 

1 

This allows moving the circumferential (Y-direction) boundaries of the inner grid. 

ENTER % PITCH FOR JBEG AND JEND (CURRENTLY=.25) 

.20 

Change the distance of the boundaries from the airfoil surface from 25% pitch to 
20%, i.e. closer to pressure and suction sides. 

ENTER IMAX (AROUND BLADE) FOR INNER GRID, FORMAT=l3 
(DEFAULTS 01, MAX=101) 

ENTER JMAX (NORM TO BLADE) FOR INNER GRID, FORMAT=l2 
(DEFAULT=21, MAX=31) 

ENTER 1 TO CHANGE DENSITY AT SURFACE 

1 

Allows changing distance from airfoil surface of 1st O-line of inner grid. 

ENTER VALUE OF DSBOD (ORIG, WITH JMX=21, IS .0005) 

.0003 



ENTER JMAX FOR TIP CLEARANCE REGION, F0RMAT=I2 
(DEFAULT=1 1, MAX=11) 


Number of O -lines in tip clearance grid. 

ENTER KMAX (RADIAL PNTS) BEFORE ADDITION, F0RMAT=I2 
(DEFAULT=25, MAX=51(NO ADDITION)) 

45 

ENTER KMAX FOR TIP CLEARANCE REGION, F0RMAT=I2 
(DEFAULT=5, MAX=15) 

13 

Optimizing K-plane distribution is a trial-and-error process; the numbers chosen 
here should be less than or equal to final number of K-planes desired. Initially , the 
program should be run through to the point at which the user has the opportunity 
to examine the distribution of the chosen number of K-planes. This should be 
repeated, with modifications made after each examination, until the distribution is 
satisfactory, at which point the program can be run through its entirety. 

CURRENT DIMENSIONS: 

NOTE THAT POINTS WILL BE ADDED IF AXIAL SPACING AT 
INLET/EXIT BOUND IS CHANGED, AND AGAIN IN OVRLAP, 

ALSO RADIAL DISTRIBUTION CAN BE MODIFIED LATER. 

AIRFOIL: 1 

OUTER: 9, 28 INNER: 1,1 RADIAL: 45 45 

AIRFOIL: 2 

OUTER: 49, 21 INNER: 101, 21 RADIAL: 45 33 
PAUSE... 

Provides an opportunity to make sure dimensions are correct before continuing. 

ENTER 1 TO INCREASE CLUSTERING ATT.E. FOR AIRFOIL 2 


This will put more grid points around the airfoil trailing edge; note this will take 
away points from the rest of the airfoil surface. 

ENTER 1 TO CHANGE ENDWALL SPACING AND/OR DIST. IN TIP 
(DEFAULT: SETA1 = SETA2 = 0.03% SPAN = 1.80000E-03) 

1 

Unless "1 " is entered, radial distance between each endwall and the next closest 
K-plane will be set to 0.03% of the span, spanwise distribution will be symmetric 
about the midspan, and the clearance region will be assumed to exist for the last 
13 K-planes, i.e., the number entered for "KMAX" for tip clearance. For most 
cases, it is better to specify the clearance region as a percentage of the total span; 
to do this, it is necessary to enter "1 " at this point. 
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CHOOSE: 1) SPECIFY TIP CLEARANCE (I.E., %SPAN) 

2) TIP CLEARANCE DETERMINED BY GRID 

1 

If option 2 is chosen, the user will be prompted for new endwall spacing, the same 
at both endwalls, and the spanwise distribution will be symmetric about the 
midspan. The tip clearance region will be defined as the top 13 K-planes. In this 
example, the span is 6 inches ; if an endwall spacing of .002 inches is chosen, the 
geometric distribution of the 45 K-planes would put the 33rd at radius=5. 82594, 
therefore the top 13 K-planes would represent 2.9% of the total span. To get the 
desired 1% clearance, option 1 must be chosen. 

ENTER PERCENT OF TOTAL SPAN FOR CLEARANCE REGION, F5.0 
(ENTRY WILL BE DIVIDED BY 100 BEFORE MULTIPLICATION) 

1 . 

This will cause the top 13 of the 45 K-planes to be disributed within the top 1% of 
the total span. 

CHOOSE: 1) K-LINES EVENLY SPACED IN CLEARANCE 

2) DENSER AT CASING THAN AT BLADE TIP 

3) DENSE AT CASING AND AT BLADE TIP 

These are the results of choosing each of the above options: 



« 

l 


2 


A 


3 


These show only the distribution in the tip clearance; the spacing through the rest 
of the span varies among the three as well. 


3 

ENTER SPACING BETWEEN CASING AND NEXT K-PLANE, F5.0 

.002 

ENTER SPACING BETWEEN BLADE HP AND NEXT K-PLANE, F5.0 

.002 

ENTER SPACING BETWEEN HUB AND NEXT K-PLANE (F5.0) 
(SPACING AT CASING IS SET TO 0.002000 

.002 
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The above are actual distances in inches. 


ENTER 1 IF YOU WANT TO ADD RADIAL LINES AT MIDSPAN 

1 

Allows filling in sparse area at midspan without taking K-planes away from 
endwalls; to determine whether this is necessary, a radial distribution should be 
generated first without excercising this option. 

ENTER NUMBER OF POINTS TO BE ADDED, 12 
(CURRENT KMX: 45) 

04 

ENTER KLO & KHI (CURRENT DIST) BETWEEN WHICH TO ADD LINES, 213 

020027 

The above will redistribute the spacing between K-planes 20 and 27 with the 
addition of four more K-planes, thus bringing the total number to 49. 

RADII OF K-PLANES ACCORDING TO ORIGINAL STRETCHING 
HAVE BEEN WRITTEN TO A FILE. 

HALT EXECUTION OR HIT ENTER TO CONTINUE 
PAUSE... 

At this point, the program should be exited and the radial distribution should be 
examined; the output file is in FORTRAN unit 30, and contains the radius defined 
for each K-plane; these can be plotted as R vs. K - here are the plots resulting from 
running this program with and without the addition of 4 lines between 20 and 27: 



without addition with addition 


Note that the distribution below K=20 is the same in both cases, as are the 
distribution above K=27 in the first plot and that above K=31 in the second. 
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AXIAL CHORD = 6.341 00 

L.E. X = 5.00000 T.E. X = 11 .341 0 

UPSTREAM BOUNDARY FROM INPUT FILE IS AT 60.8737% AXIAL CHORD 
DWNSTREAM BOUNDARY FROM INPUT FILE IS AT 37.9120% AXIAL CHORD 
(FOR VALUES LESS THAN 25, YOU MAY NEED TO REDEFINE IBEG & I END) 

ENTER 1 TO CHANGE BOUNDARIES OR RETURN 

1 

Allows changing axial boundaries from those specified in the input file. 

ENTER 1 TO SPECIFY PERCENTAGES, 2 TO SPECIFY X-VALUES 

1 

Boundaries can be specified as either percentages of axial chord, or as actual 
axial location in inches ; the latter is easier when generating grids for multi-row 
cases to ensure alignment of boundaries at the interface. 

ENTER % AX CRD AT WHICH TO DEFINE UPSTRM BOUND (F5.0) 

45 . 

ENTER % AX CRD AT WHICH TO DEFINE DWNSTRM BOUND (F5.0) 

35 . 

CURRENT AIRFOIL GRID BEING GENERATED: 2 
ENTER 0) NO CHANGE IN AXIAL SPACING 

1) CHANGE AXIAL SPACING AT INLET 

2) CHANGE AXIAL SPACING AT EXIT 

1 

For multi-row cases, increasing density of axial lines at interfaces is recommended; 
in this case, there will be a pseudo-stator upstream of the rotor, so the latter is 
considered airfoil number 2, and spacing should be reduced at its inlet. This 
involves specifying a smaller spacing than exists when the domain is divided evenly 
by the number of axial points in the grid, and adding an appropriate number of 
axial points to compensate for the reduced spacing at one end. 

ENTER AXIAL SPACING AT BOUNDARY (FI 0.0) 

NOTE: EVEN SPACING IS 0.237787 
.065 

ENTER NUMBER OF NEW POINTS TO BE ADDED (il) 

(Suggestion: 4 

4 

Determining the proper spacing and number of additional points may require some 
trial and error, but one-third to one-quarter of the even spacing is a good start; 
note that the total number of axial points must not exceed program dimensions, 
also, inlet spacing for an airfoil must match exit spacing of the preceding one. 

ENTER 1 TO FURTHER CHANGE SPACING, OR RETURN 
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ENTER VALUE OF YFRAC(I), OR RETURN FOR -.5 

This refers to the distance between the airfoil camber line and the lower 
circumferential boundary of the outer H-grid in terms of the fraction of the pitch. 
In almost all cases, the default (i.e. -.5 or 50% pitch below the camber line) will 
work. 

NOTE: IB,IE,JB,JE: 15 46 7 15 

CHECK TO MAKE SURE THESE ARE OK BEFORE INITIALIZING. 

PAUSE... 

At this point, these values should be recorded, and upon completion of the grid, 
check to make sure the outer H-grid lines defined by I— IB and I— IE for 
JB-1 <J>JE+1 , and by J=JB and J=JE for IB-1 <I>IE+I , are contained entirely 
within the boundaries of the inner O-grid; if not, the program should be re-run 
with appropriate adjustments made to IBEG & IEND andlor JBEG & JEND. 
ENTER 1 TO SKIP ELLIPTIC GRID GENERATOR 
1 

Recommend skipping elliptic grid generator, which requires a significant amount 
of CPU, until the final run of the program. 

ENTER K-INDEX TO BE USED TO DETERMINE PITCH, 12 (KMX: 49 

25 

This is required only for a 3D planar case; the pitch will be the same for all 
K-planes, and will be equal to 2* PI* (radius at chosen K)!(number of blades). 

In this case, K—25 is approximately at midspan, at which the radius is 27 inches, 
therefore the pitch will be 2*PI*27I28 = 6.059 inches. 

ENTER PERCENT 2ND ROW BX UPSTREAM OF 2ND ROW AT WHICH 
PSEUDO-AIRFOIL INLET WILL BE DEFINED 
(NOTE: ROTOR INLET BOUNDARY IS DEFINED AT 45.0000 

60 . 

The pseudo-stator inlet will be 15% rotor bx upstream of the rotor inlet. 

ENTER FACTOR OF 2ND ROW GAP TO DEFINE PSEUDO-AIRFOIL GAP 
(E.G., FOR A 3-VANE, 4-BLADE CASE, ENTER 1.333) 

1.333 

The pseudo-stator grid pitch will be 4/3 that of the rotor . 

ENTER 1 TO SHIFT PSEUDO- GRID, RETURN FOR ALIGNMENT W/2ND ROW 

1 

For viewing purposes, most useful for 1/1 ratios; shifts pseudo-airfoil grid in 
circumferential direction. 

ENTER % PITCH TO SHIFT GRID (F5.0) 

(CAN BE NEGATIVE) 
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- 10 . 

WRITING PLOT3D FILE ... 

% 


The above run produces the following grid: 



View in radial direction View in axial direction 

Note the inner O-grid is not finalized. At this point, if all other aspects of the grid are 
satisfactory, the program should be run again without skipping the elliptic grid 
generator. 

In summary, the steps to generate a grid are as follows: 

1. Create the input file containing the airfoil(s) geometry in unit 2. 

2. Invoke the program and, for a 3D case, exit after radial distribution 
has been written out. 

2a. (3D only) Repeat step 2 until desired radial distribution is achieved. 

3. Run program through, excercising option to skip elliptic generator. 

4. Repeat step 3 until desired inner/outer grid relationship is achieved. 

5. Run program through, do not skip elliptic grid generator. 

The following is a list of the subroutines and their functions. 
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MAIN: 


prompts user for general information about the case; calls: 


DIMEN: 

GNRLAF: 

GEOM: 

RSPACE: 

SPLINE: 

CSPLDM: 

CAMBER: 

TSPLIN: 

ALGEB: 

OBOUND: 

ELLIP: 

TRI: 

REDIST: 

OVRLAP: 

PSEUDO: 

ROTLID: 

UNWRAP: 

WRPOUT: 

GEOMOT: 

PLOT3D: 


prompts user for grid dimensions 

main grid generation routine, called for each airfoil; calls: 
reads in airfoil geometry from input file 
for 3D cases, prompts user for spanwise grid distribution 
fits a spline to the input airfoil coordinates, defines distribution of 
grid lines normal to airfoil surface based on curvature 
calculates derivatives of X and Y wrt S 
generates airfoil camber line in X-Y coordinates 
calculates curvature of outer grid boundaries 
generates outer H-grid 
generates outer boundary of the inner O-grid 
elliptic grid generator, generates the inner O-grid 
solves the tridiagonal matrix calculated in ELLIP 
redistributes inner grid O-lines based on chosen surface spacing 
for multi-row cases, generates patched boundary at each row 
interface 

for "pseudo-airfoil" cases, generates a small grid upstream of first 
airfoil grid, representing the exit of the previous airfoil 
for cases with tip clearance, generates tip clearance grid 
for planar cases, "unwraps" the annular grid about the X-axis 
for radial cases, "unwraps" the outer H-grid about the Z-axis 
writes grid and case information to file for input to initializer 
writes grid to PLOT3D-readable file 


Miscellaneous routines: 


SPACE 1, 

SPACE2: geometric point distributors 

PUTXYZ, 

GETXYZ: puts and gets coordinates from temporary storage files 

OPNCLO: opens and closes temporary storage files 

The format of the input file is described here. All of these lines must be present 
for every full airfoil for which a grid will be generated. ( No input is required for a 
pseudo-airfoil, its grid will be determined by grid with which it interfaces.) 
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Card# 

Columns 

Format 

Variable 

Description 

1 

1-80 

A80 

(dummy) 

Title card to separate the airfoils; 
not used by the program, but line 
must be present. 

2 

1-5 

15 

NBLADE 

Total number of blades in this 
airfoil row; wll be used in 
conjunction with radius to 

Card# 

Columns 

Format 

Variable 

Description 





determine periodic boundaries 
of solution domain. 

3 

1-20 

G20.8 

XU 

Axial location of upstream bound. 

3 21-40 

(Note that XU and 

G20.8 XD Axial location of downstrm bound. 

XD maybe be changed during program execution.) 

4 

1-5 

15 

NSPANS 

Number of spanwise locations 
(radii) present to define airfoil 
geometry; MIN = 2, MAX = 10 

4 

(for N= 

6-10 

=1. NSPANS:) 

15 

NSTR 

Number of X-Ycoordinates given 
to define airfoil geometry at each 
spanwise location; MAX = 251 

5 

(for 1= 

1-20 

l.NSTR:) 

G20.8 

R(N) 

Radius at which the following 
NSTR lines define the airfoil 
geometry 

6 

1-20 

G20.8 

X(N,I) 

X-coordinate of Nth point at Ith 
spanwise location 

6 

21-41 

G20.8 

Y(N,I) 

Y-coordinate of Nth point at Ith 
spanwise location 


Card type #6 is then repeated to make a total of NSTR lines containing X-Y values. 
Note the X-Y coordinates must proceed clockwise around the blade 


The Card#5/Card#6 combination is repeated to make a total of NSPANS groups of 
NSTR+1 lines, i.e., one line formatted as Card #5 followed by NSTR lines formatted 
as Card#6. Each group defines the airfoil geometry at the specified radius. The first 
radius must be at the hub, and the last must be at the casing. Note that the program 
will linearly interpolate the coordinates at the given radii to define the geometry at 
25%, 50%, and 75% span, and it is those, in addition to the hub and casing geometry, 
that will be used to generate the grid. 
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a 


The entire format (i.e., beginning with a title card) is repeated for subsequent airfoils. 
The following is the input file accompanying the preceding sample run: 


LSRR rotor 
28 

1.14 
5 136 

0 . 24000000E+02 
0 . 500000 00E+01 
0 . 500132 94E+01 
0.50052996E+01 
0.50118799E+01 
0 . 50210295E+01 
0 . 50326691E+01 


13.745 


0.28807001E+01 
0 . 28503103E+01 
0 . 28201504E+01 
0 . 27904596E+01 
0 . 27614498E+01 
0 . 273334 98E+01 

of this format) 


0.33143501E+01 
0 . 32839804E+01 
0 . 32538404E+01 
0 . 32241602E+01 
0 . 31951704E+01 
0 . 31670904E+01 

of this format) 


0 . 37861996E+01 
0 . 37558098E+01 
0 . 372564 98E+01 
0 . 36959400E+01 
0 . 36669302E+01 
0 . 36388197E+01 

of this format) 


0 . 42581501E+01 
0.42277298E+01 
0. 4197529 8E+01 
0 . 41677999E+01 
0 . 41387596E+01 


(130 more lines 

0 . 25500000E+02 
0.50000000E+01 
0.50013294E+01 
0.50052891E+01 
0.50118694E+01 
0.50210094E+01 
0.50326500E+01 

(130 more lines 

0.27000000E+02 
0 . 50000000E+01 
0 . 50013294E+01 
0 . 50052 996E+01 
0 . 50118799E+01 
0 . 50210295E+01 
0 . 5032 6691E+01 

(130 more lines 

0 . 28500000E+02 
0 . 50000000E+01 
0 . 50013294E+01 
0 . 50052 996E+01 
0 . 50118895E+01 
0 . 50210495E+01 
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0 . 50327091E+01 

(130 more lines 

0 . 30000000E+02 
0.50000000E+01 
0.50013294E+01 
0 . 50 052 99 6E+ 01 
0.50118895E+01 
0 . 50210400E+01 
0.50326796E+01 

(130 more lines 


0 . 41106195E+01 
of this format) 


0 . 46918201E+01 
0 . 46614199E+01 
0 . 46312504E+01 
0 . 46015396E+01 
0 . 45725203E+01 
0 . 45444098E+01 

of this format) 
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Flow Initializer 


The purpose of the flow initializer is to define an initial guess at the solution of 
the flow to give the flow solver something to start with. In addition, it is in this 
program that the boundary conditions that will be applied in the solver are defined. It 
requires, as input, the aerodynamic conditions at the inlet and exit of the solution 
domain; also, the number of airfoils in each row for which a solution is desired, the 
speed of rotation, Reynolds number, and Prandtl numbers for the case. The 
computational grid, i.e. the file "GEOM", output from the GRIDGEN, is the only 
other input, and there is no interactive user input. Output is in two forms: l)a binary 
"restart" file to be input to the flow solver, 2) a PLOT3D file for viewing the initial 
flowfield and boundary conditions. 

The method of initialization is slightly modified from Rai’s original version. 
Static pressure and density are linearly interpolated between the inlet and exit of the 
solution domain, axial and radial velocities are constant across the domain, and 
tangential velocities are linearly interpolated between the inlet and exit of each airfoil 
based on the general curvature of the airfoils. A no-slip condition is then applied to 
all surfaces, which was not done by Rai. This is done by defining the (relative) 
velocity on each surface to be zero, then linearly interpolating to a pre-defined 
distance from the surface along grid lines that are normal to the surface. Another 
added feature not present in Rai’s original code is the dependency of boundary values 
on radial and tangential location. This allows specification of an incoming endwall 
boundary layer (if desired) in three dimensions, and/or an incoming wake to simulate 
the presence of an airfoil upstream of the inlet to the solution domain, in either two or 
three dimensions. The same flexibility applies to the exit plane, where measured 
static pressure data may be held across the plane if the case so warrants. 

The following is a list of subroutines and their functions: 


MAIN: 

GEOMIN: 

AEROIN: 

BCDEF: 

INIT: 

INITOT: 

PLOT3D: 


main driver, calls: 

reads file "GEOM", containing grid coordinates and case info, 
reads file containing aerodynamic and viscous parameters 
interpolates and nondimensionalizes inlet and exit flow properties 
interpolates between inlet and exit to initialize flowfield, imposes 
no-slip condition on surfaces 

writes file "INIT", restart file for input to flow solver RAI3DS 
writes grid and initial flow in PLOT3D format for viewing 


Miscellaneous routines: 
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SHIFT: 

PUTXYZ, 

GETXYZ: 

PUTQ, 

GETQ: 

OPNCLO: 


applies periodic shift to airfoil grids for multi-blade cases 

puts and gets grid coordinates from temporary storage files 

puts and gets flow properties from temporary storage files 
opens and closes temporary storage files 


The format of the aerodynamic input file is described here. There are two 
possibilities for defining the inlet and exit conditions. The first is to specify inlet 
aerodynamic conditions and exit static pressures at specific radii. These will then be 
linearly interpolated onto the radii of the computational grid, and will not vary in the 
circumferential (tangential) direction . How properties will then be calculated from 
the aerodynamic conditions and nondimensionalized according to the scheme 
expected by the solver. The other option is to input dimensional values of the flow 
properties themselves (density, velocity components, and pressure) at every 
computational grid point in the inlet plane, along with dimensional static pressure at 
every grid point in the exit plane. These will then be nondimensionalized according 
the the scheme expected by the solver. The latter option will allow specification of an 
incoming wake, but it is much more difficult, particularly in three dimensions. It 
requires interpolation of the flow properties from some other source, whether it be 
experimental data or output from a flow solver, etc., onto the grid coordinates. A 
method of accomplishing this interpolation is not provided in the RAI3DS system. 
This file must be in FORTRAN unit 2; here is its format: 


Card# 

Columns 

Format 

Variable 

Description 

1 

1-80 

A80 

(dummy) 

Title card to identify case; not used 
by program but must be present. 

2 

1 

11 

IV ARY 

=0: aerodynamic input is specified 
at radii, does not vary across pitch; 
=1: flow properties are defined at 
grid points. 

3 

1-2 

12 

NBLDS(l) 

Number of blades in row 1 for 
which solver will calculate flow 
(MAX = 4). 

3 

3-4 

12 

NBLDS(2) 

Number of blades in row 2 for 
which solver will calculate flow 
(MAX = 4). 
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Card# 


Columns Format Variable Description 


3 5-6 12 NBLDS(3) Number of blades in row 3 for 

which solver will calculate flow 
(MAX = 4). 

(Note: NBLDS(1 )+NBLDS(2)+NBLDS(3) must not exceed 7.) 


4 

1-10 

FI 0.0 

RPM 

Speed of rotation in revolutions 
per minute (=0.0 for stationary 
cases); can be negative. 

4 

11-20 

FI 0.0 

GAMMA 

Specific heat ratio. 

4 

21-30 

FI 0.0 

CP 

Specific Ipit at constant pressure 
in (ft /sec /°R). 

4 

31-40 

FI 0.0 

REYNIN 

Reynolds number per inch based 
on inlet conditions. 

4 

41-50 

F10.0 

PRKIN 

Laminar Prandtl number. 

4 

51-60 

F10.0 

PRTUR 

Turbulent Prandtl number. 

If IVARY= 

=0: 




5 

1-2 

12 

KUP 

Number of radii at which inlet 
aerodynamic conditions will be 
defined; MIN=1, MAX=50. 

6 

1-10 

FI 0.0 

RADUP 

Radius at which the conditions on 
this line apply. 

6 

11-20 

F10.0 

P0INF 

Absolute total pressure at inlet, 
in (lb/in ). 

6 

21-30 

FI 0.0 

T0INF 

Inlet total temperature in (°R). 

6 

31-40 

FI 0.0 

AMACH 

Absolute Mach number at inlet. 

6 

41-50 

F10.0 

ALPHA 

Absolute inlet flow angle in 
degrees,defmed as tan’ (V/U), 
where V=tangential velocity and 
U=axial velocity. 

6 

51-60 

F10.0 

PHI 

Absolute inlet flow angle in 
degrees, defined as tan (W/U), 
where W=radial velocity and 
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Card# 


Columns Format 


Variable Description 


U=axial velocity. 

Card 6 is then repeated to make a total of KUP lines. 


1-2 

12 

KDW 

Number of radii at which exit 




static pressure will be defined; 
MIN=1, MAX=50. 

1-10 

FI 0.0 

RADDWN 

Radius at which the static pressure 
on this line applies. 

11-20 

FI 0.0 

PDWN 

Exit static pressure in (lb/in 2 ). 


Card 8 is then repeated to make a total of KDW lines. 


If IVARY=1: 

Forn=lJSfBLDS(l): 

For k=l,KMX, where KMX=number of radial planes in grid: 

For i=l JMX(lk where JMX(l)=number of tangential points in 1st row grid: 


5 

1-15 

G15.8 

RINF(n j,k) Inlet density in (lb-sec 2 /in 4 ). 

5 

16-30 

G15.8 

UINF(nj,k) Inlet axial velocity in (in/sec). 

5 

31-45 

G15.8 

VINF(n j ,k) Inlet tangential velocity in (in/sec); 

note this is tangential velocity for 
annular cases, Y-vel. for planar. 

5 

46-60 

G15.8 

WINF(nj,k) Inlet radial velocity in (in/sec); 

note this is radial velocity for 

annular cases, Z-vel. for planar. 

5 

61-75 

G15.8 

PINF(n,j,k) Inlet static pressure in (lb/in 2 ). 


(This card will appear a total of NBLDS*KMX*JMX times) 


For n=l,NBLDS(last row): 

For k=l,KMX, where KMX=number of radial planes in grid: 

For i=lJMX(last row), where JMX=number of tangential points in last row grid: 

6 1-15 G15.8 PDWN(n,j,k) Exit static pressure in 

(lb/in). 

(This card will appear a total of NBLDS*KMX*JMX times) 
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Note that if the case is planar, either two- or three-dimensional, the geometry 
will have been translated in the grid generator so that the first K-plane has a constant 
Z-value of 0.0. The values of RADUP and RADDWN, however, must correspond 
with the original radii input to the grid generator. The following is a sample input file 
corresponding to the sample run of the grid generator previously shown: 

LSRR rotor 
0 

1 1 


0.00 

1.4 

6007.120 

41000.0 

0.72 

0.90 

\ 

24.00000 

14.60446 

540.00000 

0.15500 

68.2272 

0.20240 

24.30000 

14.68998 

540.00000 

0.17721 

67.32001 

0.44000 

24.40021 

14.68890 

540.00000 

0.18047 

66.87000 

- 0.89000 

24.49980 

14.68432 

540.00000 

0.18016 

66.62000 

- 1.21000 

24.70020 

14.67641 

540.00000 

0.17641 

65.99001 

- 1.55000 

24.89999 

14 . 67336 

540.00000 

0.17474 

65.78999 

- 2.01000 

25.30020 

14 . 69601 

540.00000 

0.18069 

67.33000 

- 2.65000 

25.80000 

14.69550 

540.00000 

0.17372 

67.05000 

- 0.26000 

26.39999 

14.69499 

540.00000 

0.17143 

67.39000 

- 0.77000 

27.00000 

14.69397 

540.00000 

0.16894 

68.39000 

- 2.95000 

27.60120 

14 . 69056 

540.00000 

0.17413 

68.27000 

- 3.12000 

28.20000 

14.68664 

540.00000 

0.16231 

67.27000 

- 3.44000 

28.50121 

14.68991 

540.00000 

0.16094 

67.13000 

- 1.47000 

28.70100 

14.69310 

540.00000 

0.16065 

67.39000 

- 1.37000 

28.90021 

14 . 69557 

540.00000 

0.16098 

67.22000 

- 1.15000 

29.10120 

14.69477 

540.00000 

0.16118 

67.14000 

- 1.07000 

29.30099 

14.69441 

540.00000 

0.16004 

68.41000 

- 0.70000 

29.49120 

14.69419 

540.00000 

0.15836 

71.52000 

- 0.24000 

29.71730 

14.63597 

540.00000 

0.14163 

72.33858 

- 0.13300 

30.00000 

14 . 60228 

540.00000 

0.12944 

73.58221 

- 0.07090 


2 

24.00000 14.2018 

30.00000 14.21483 


Note that the boundary values and initial flowfield should be in absolute frame 
and Cartesian coordinates. This should be checked using PLOT3D before running the 
flow solver. When viewing the flowfield, also note that density has been nondimension- 
alized by the gap-averaged inlet density at midspan, call it RNDF; pressure, by the 
gap-averaged inlet pressure at midspan, call it PNDF; velocity components have been 
nondimensionalized by V(PNDF/RNDF). 
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Row Solver 


The flow solver’s main input is a "restart file", which contains all geometric and 
aerodynamic information about the case, including flow properties at each node that 
were calculated at the end of the solver’s previous run. In addition, a small amount of 
input affecting a particular run, such as number of time steps, size of time step, 
frequency of printing convergence information, etc., is required. (For two-dimensional 
cases, a file defining h-ratio is also necessary.) For each time step, the solver runs 
through the solution algorithm for each airfoil grid separately, interpolating at grid 
boundaries, explicitly enforces boundaiy conditions, and records convergence 
information. For unsteady cases, this is repeated (usually twice) before the rotating 
grid(s) is moved relative to the stationary grid according to the size of the time step. 
After the last time step, the restart file is updated with the latest flowfield and output. 
For unsteady cases, there is additional output consisting of files containing pressure 
distributions on airfoil surfaces, and total pressure across the exit plane, which have 
been written at specified time intervals. This allows analysis of unsteady behavior as 
a function of time, in addition to the instantaneous solution recorded in the restart file. 

The program solves the unsteady, three-dimensional, compressible Navier-Stokes 
equations, supplemented by an equation of state and an energy equation, all of which 
are first transformed from Cartesian coordinates into a (£, T|, q) coordinate system. In 
this coordinate system, t, is the direction tangent to the airfoil surface, rj is the 
direction normal to the airfoil surface, and <; is the direction normal to the hub surface. 
The thin-layer assumption is made, that is, the viscous terms evaluated as derivatives 
in the direction tangent to the body surface are assumed to be negligible. The 
transformed equations then take the following form: 

Q + E+ F +G = Re *(S + T) (1) 

, * % <; ? 
where 

Q = Q/J 

m,%) = G,Q + % E + IF + %G)IJ 

F(Q,i\) = (t\Q + t^E + t) F + t\G)IJ (2) 

G( Q,$) = (g<2 + sj E + <;$ + zG)/J 

Q is the dependent variable vector, (p, p u, pv, pu\ e) in Cartesian coordinates, where 
p is density, u, v, and w are velocities in the x-, y-, and z-directions, respectively, and e 
is energy; E, F, and G are flux vectors in the x-, y-, and z-directions, respectively; S 
and T are the viscous flux vectors in the T|- and ^-directions, respectively; J is the 
Jacobian of the coordinate transformation; Re is Reynolds number. The integration 
scheme used to solve the set of equations is described in detail in Ref. 1 . It’s final 
form is as follows: 
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where 

A = (dEldQ) ± 
/ = ( dF/dQ ) ± 
C = (dG/dQ) 
M = (dS/dQ) 
N = (dT/dQ) 


(3) 


(4) 


and A, V, and 8 are forward, backward, and central difference operators, respectively. 


The E F ,G ,S , and T are numerical fluxes consistent with 

i+l/2J,k’ ij+mjk ij,k+H2 ij+l/2,k ijMl/2 . , /D r a \ 

their corresponding physical fluxes, and are evaluated using Roe s scheme (Ret. 4). 


The superscript n refers to the time step, and the superscript p refers to the iteration 
per time step. The subscripts ij,k refer to the grid node. As mentioned before, the 
derivation of this scheme and it’s limitations are fully described in Ref. 1; it is 
presented here only for reference. 

The following is a general list of subroutines and their functions; a detailed 
version can be found in Appendix 1. Key features and modifications made to some 
of these routines, i.e., enhancements to Rai’s original ROTOR3 code, are documented 
in some detail following the list. They are categorized according to the general 
feature which was modified, which may have required implementation in several 
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routines, as opposed to simple listing of the changes made to each routine. 

RAI3DS (main routine): reads indicators to dictate run: calls: 

RESTIN: reads initial flow file or restart file from previous RAI3DS run 

JACBIN : calculates Jacobians for the variable transformation 

TRIANG: calculates interpolation factors at inner/outer interface 

HDATA: (2D only) reads H-ratio input data 

EIGEN: (steady only) calculates "time" step 

(for each iteration) . 

CONTRL: controlling routine for each solution iteration; calls: 

GMOVE: (multirow only) moves rotating grids relative to stationary grids 

PGRID: defines solution above and below periodic bounds for use in RHS 

RHS: control routine for computing the right hand side of Eq. (3); calls: 

SRINT: (outer grids only) interpolates at inner/outer interface 

MUKN: calculates kinematic viscosity (Sutherland’s law) 

with modifications) 

, and Cj in Eq. (3) 

P in Eq. (3) 

STCONT: (2D only) adds stream tube contraction terms in place of Cj 
LHS: control routine ^or computing the left hand side of Eq. (3) and 

solving for Q ; calls: 

SMATRX:computes the matrices A ,B , and C in Eq. (3) 

VMAT: computes the matrices M and N in Eq. (3) 

BTRI: inverts the matrix equation (3), solving for Q 

CORREC: imposes explicit boundary conditions at all boundaries 

CONVRG: calculates change in solution from previous iteration . 

WRTPS: writes out pressure distribution at midspan 

WRTEXT : writes out exit total pressure across pitch at midspan 

RESTOT: writes out restart file containing updated solution 

2D or 3D solution : RAI3DS has the capability for either two- or 
three-dimensional solutions. Its algorithms work with three dimensions, and for this 
reason, even if a two-dimensional solution is desired, a three-dimensional grid is 
required. The grid for a two-dimensional solution must consist of three identical grid 
planes at three separate Z-locations. A maximum K-index of 3 will signal the solver 
to compute the solution at only the middle plane and copy it to the planes above and 
below it. This will ensure zero flux in the Z-direction. In addition, the solver wjll 
skip evaluation of terms with respect to Z, i.e., vectors G and T and matrices C and N 
in Eq (3), saving on computation time. (Affected subroutines: CORREC, all RHS 


MUTR: calculates eddy viscosity (B aid win-Lomax 

FLUXR: computes the numerical flux vectors E?, F* 
VFLUX: computes the viscous flux vectors and 7 
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and LHS; also, routine MUTR2D will be called instead of other MUTR routines.) 

Stream tube contraction terms : Some three-dimensional effects have been 
shown to be effectively captured in a two-dimensional solution by adding stream tube 
contraction terms (Rangwalla, Ref. 5). The calculation of these terms and their 
incorporation in the solution algorithm have been added to RAI3DS. (Affected 
subroutines: all RHS routines and VFLUX, plus addition of HDATA and STCONT.) 

Planar or annular configuration : Although RAI3DS always works in Cartesian 
as opposed to cylindrical coordinates, the configuration of the geometry may be either 
a plane or an annulus. Obviously, a two-dimensional case must be planar, but a 
three-dimensional case may have either configuration. (Affected subroutines: 
GMOVE, PGRID, CORREC; logical variable PLANE is indicator.) 

Symmetry condition : For axisymmetric three-dimensional cases, more grid 
resolution can be obtained by using the solver’s full grid capacity for the region 
between the hub and the midspan of the airfoil and imposing a symmetry condition at 
the uppermost K-plane. This, in effect, doubles the spanwise resolution. Obviously, 
this is applicable only to planar three-dimensional cases with no tip clearance. It has 
been demonstrated with such a case, however (Ref. 6), that finer spanwise grid 
resolution significantly improves the solver’s capability for predicting 
three-dimensional effects. For details on this feature and instructions on how to use it, 
see Appendix 2. (Affected subroutines: all RHS, all LHS, and all MUTR routines, 
CORREC; logical variables HAFSYM and FULSYM are indicators.) 

Multi-blade, multi-row capability : Current limitations in the solver are four 
blades per row and three total rows. (Affected subroutines: GMOVE, PGRID, outer 
LHS routines, CORREC; variables NBLADS and NROWS are indicators.) 

Radial flow capability : While ROTOR3 was applicable only to axial flow 
machines, in which rotation occurs in the plane normal to the direction of the 
incoming flow, RAI3DS can be used to solve cases in which the rotation occurs in the 
same plane as the incoming flow vector. Specifically, for axial flow cases, rotation is 
about the X-axis, and for radial cases, rotation is about the Z-axis; in both cases, flow 
convects along the X-axis. The current version of the code can be used for impellers 
without tip clearance or inlet guide vanes/diffusers. (Affected subroutines: 

CORREC, GMOVE, PGRID, TRIANG, SHIFT, EIGEN, all RHS, all LHS, all 
MUTR; logical variable RADIAL is indicator.) 

Inlet/exit boundary conditions : The initial values of entropy, pitch angle, yaw 
angle, and total enthalpy are held at the inlet, accompanied by a Reimann invariant 
that is extrapolated from the interior. Holding both the enthalpy and entropy results in 
holding total pressure, which is a property commonly specified by design engineers at 
the inlet. In addition to allowing total pressure specification, this is a less stiff 
condition than that in ROTOR3, which held initial entropy, tangential and radial 
velocities, an initial Reimann invariant, and a Reimann invariant extrapolated from the 
interior. The relaxed conditions reduce the occurence of pressure waves, which may 
slow the convergence of the calculation in some cases. Specified static pressure is 
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held at every point in the exit plane, modified from the radial equilibrium condition in 
ROTOR3. Both the inlet and exit conditions may vary in both the tangential and 
radial directions, allowing specification of more general conditions than were possible 
in ROTOR3. The derivation of the inlet boundary condition is presented in Appendix 
3. (Affected subroutine: CORREC.) 

Surface boundary conditions : Capability of specifying wall temperature or heat 
flux on solid boundaries was added to ROTOR3. These options may be chosen, or an 
adiabatic wall condition may be used. Hie implementation of specified wall 
temperature was derived by Griffin (Ref. 7), and the implementation of specified heat 
flux is presented in Appendix 3. (Affected subroutines: LHS, CORREC, plus 
addition of BCSURF; ISURF is indicator.) 

Turbulence model : The Baldwin-Lomax turbulence model in ROTOR3 was 
modified to account for surface curvature, three-dimensionality and Coriolis force. 
These modifications, deduced by interrogating transport equations for the turbulent 
kinetic energy and total Reynolds shear stress, are similar to those termed by various 
investigators as ’algebraic stress models’ (Ref. 8). A transition model has also been 
implemented that accounts for the influence of upstream airfoil wakes on the onset 
and the extent of the transition on the downstream airfoils. This model is primarily 
applicable on the suction sides of airfoils operating at moderate Reynolds number. 
(Affected subroutines: all MUTR routines.) 

Pseudo-airfoil : This capability was added to allow solution for a rotor moving 
past a stationary wake without having to solve for the upstream stator. Conditions are 
specified at the inlet of a small H-grid, which remains stationary and is treated as an 
"outer" grid by the code. The following rotor moves past the H-grid as it would a 
stator grid. (Affected subroutines: CONTRL, CORREC, JACBIN, TRIANG, 

RHSSO; variable IPSEUD is indicator.) 

To create the necessary input and execute the flow solver, the user should run 
the shell script, RAUOB. This script is run on die Iris and will prompt the user for all 
information necessary for defining the input files, setting up the file environment, and 
executing the flow solver. The output from the script is a UNICOS jobstream, which 
is then submitted to the Cray. Presented below is a sample run of the shell script 
corresponding to the three-dimensional pseudo-stator/rotor case for which a sample 
run of the grid generator was previously shown. Note that in order for the script to 
execute, the file "skel.job" must exist in the current directory. 
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% RAI JOB 


********************************************************************** 


* RAI3DS 

* * 
********************************************************************** 


This shell script creates a UNICOS jobstream for running RAI3DS on the NASA/MSFC 
Cray. 

You will be asked questions about your job, and the file "skel.job" will be edited. 

The file "skel.job" must be in the current directory. 

(Enter "q" to quit or Carriage Return to continue.) 

Please enter Cray userid 

ckabl96 

Please enter Cray password 

passwd 

Please enter job name (7 characters max) 

lsrr 

Information about the case size is required to request SSD space 

Please enter number of airfoil rows in the case 

2 

Note this includes the pseudo-stator. 

2D OR 3D? (Enter 2 OR 3) 

3 

How many airfoils in row 1 ? 

1 

Tip clearance in row 1? (y or n) 

n 

How many airfoils in row 2 ? 

1 

Tip clearance in row 2? (y or n) 

y 

Please enter KMAX 

49 
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Please enter numer of k-planes in tip clearance 

13 

Please enter name of input restart file residing in $HOME on the Cray 

lsrrinit 

The script is set up to copy the file ’ $HOMEI lsrrinit’ into FORTRAN unit 8, 
which the solver expects to contain the input restart file. This can be changed 
by editing the UNICOS jobstream if the input file does not reside in the user’s 
home directory on the Cray. 

Please enter name of output file to be written to $HOME, or 
Carriage Return if you do not want to save the file on the Cray. 

(Either way, the file will be disposed back to the IBM.) 

Isrr200 

Please enter number of steps per rotation through one rotor passage. 

3000 

This determines the time step. The number is case-dependent, and for some 
cases in which the grid is particularly dense, the solver may blow up if the time 
step is too large. It is safest to use at least 2000 steps per rotor passage. For 
steady-state cases, the time step is based on grid density and the flowfield, and 
can be increased to speed up convergence of the solver. 

Please enter number of iterations per time step 
3 

This is usually 3 for unsteady cases and 1 for steady-state cases. 

Do you want to change frequency of eddy viscosity calculation? (y or n) 

y 

If yes, eddy viscosity will be calculated at the first time step, then only as often 
as is specified. Between updates, the most recently calculated value is used at 
each grid point. This saves on CPU, since eddy viscosity calculation is very 
time-consuming, and does not adversely effect the solution as long as the 
viscosity is updated reasonably often (every 50 or 100 time steps). 

Please enter iteration interval at which eddy viscosity will be calculated 
(e.g., for "50", eddy viscosity will be updated every 50 iter) 

100 

Please enter iteration number at which to impose this frequency 
(e.g., for "500", eddy viscosity will be calculated at every iteration 
until number 500, after that, every 100 iterations 

1000 

Do you want to force transition? (y or n) 

y 

This requires knowing where the transition points are on the airfoil in terms of 
grid indices. The turbulence model will calculate the natural transition point 
using empirical correlations, which may be inaccurate. If the transition point is 
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known (from experimental data), best results are obtained by specifying it. 

Please enter l-indices of inner O-grid between which eddy viscosity 
will be turned on (same for all rows) 

(e.g., for "37 81 ", eddy viscosity will be set to zero for 1=1 ,37 and for 1=81 ,IM AX) 

56 92 

Please enter K-planes between which calculation will be transitional 
(e.g., for "7 25", calculation will be fully turbulent 
around entire blade for K=1,7 and for K=25,KMAX). 

Hit Carriage Return for full span transitional. 

Enter: 

0) Surface boundary condition will be Ist-order, adiabatic 

1) Surface b.c. will be Ist-order, specified wall temp 

2) Surface b.c. will be 2nd-order, specified wall temp 

3) Surface b.c. will be 2nd-order, specified heat flux 

3 

See Appendix 3 for descriptions of these. 

Enter heat flux on blade, hub, casing, separated with spaces; must be reals with 
not less than 1 digit before the decimal point and 5 digits after the decimal point 
(e.g., 0.41500); 

see Users Manual for instructions on how to calculate these in proper units. 

0.41500 0.52500 0.52500 
See Appendix 3. 

Enter: 

0) No line of symmetry will be imposed 

1) Restart file contains full span 

(program will solve for K=1 to midspan, then impose symmetry) 

2) Restart file contains half span 

(program will solve for K=1 to KMAX with modified b.c. at KMAX) 

0 

See Appendix 2 for details. 

How many time steps? 

200 

How often do you want convergence info written out? (every _ time steps)? 

10 

This consists of the maximum change in energy over each grid for each airfoil 
and its grid indices. It should be noted that this information can be used to 
insight into the convergence rate of a case and/or to locate trouble spots in the 
flowfield. To decide whether a case is converged, however, requires thorough 
examination of the flowfield, particularly surface static pressures and exit total 
pressure contours. For this reason, there is no predefined convergence 
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criterion built in to the solver. 

How often do you want midspan Ps dist. written out? (0 = not at all)? 

50 

This is useful in unsteady cases for plotting pressure amplitudes. The 
(nondimensional) surface static pressures at midspan are written as a function 
of time step, I -index of the O-grid, and percent axial chord downstream of the 
leading edge. This format can be modified in the routine WRTPS to produce 
proper input formats for various plotting packages. In this case, the 
information will be written out every 50 time steps. 

How often do you want exit Pt written out? (0 = not at all)? 

50 ! 

Same as above, except total pressure across the exit plane at midspan will be 
written out; format is in routine WRTEXT. 

Please enter CPU time in seconds 

5000 

This is best determined by trial and error; a typical 2D cascade requires less 
than 112 second per iteration, whereas a 3D stage may require over 30 seconds 
per solution iteration (90 seconds per time step), depending on grid size. 

The file Isrr.job has been created. 

% 


The above run produces the following jobstream: 

# USER=ckab196 PW=passwd 

# QSUB -r Isrr 

# QSUB -eo 

# QSUB -Im 4.0mw 

# QSUB -It 4990 

# QSUB -IT 5000 

# 

# CREATE TEMPORARY DIRECTORY AND TEMPORARY SSD 

# 

set -x 

cd $TMPDIR 
SSD=*tmpdir /ssd‘ 

# 

# CLEAN-UP IF JOB TERMINATES ABNORMALLY 

# 
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trap "cd ..; exit" 0 2 3 15 26 
# 

# BUILD SEGLDR INPUT COMMAND FILE AND LINK 

# fetch the object file from the Cray; to create the object file, 

# compile the source using cft77 with the -b option 

# 

fetch slvr.o -f TR -t , dsn=CKAB196.UNIC.OBJECT.RAI3DS,DISP=SHR’ 
# 

cat >seg.input«eofseg 

bin=slvr.o 

abs=a.slvr 

xfer=RAI3DS 

eofseg 

segldr seg.input 
# 

# DEFINE FILE 3 - main input file 

# 

cat >fort.3«eof3 
25 

002000001 00005000050 

2000300003 

01.0099999 

01001000 

056092 

001049 

3 

0000.41 5000000.525000000.52500 
0 

eof3 

# 

# DEFINE FILE 4 - hratio values (this file is not used for 3d cases) 

# 

cat >fort.4«eof4 


nvals: 

2 

X 

h 

-100. 

1.0 

100. 

1.0 

nvals: 

2 

X 

h 

i 

o 

p 

1.0 

100. 

1.0 

eof4 


# 
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# DEFINE FILE 8 (input restart file) 

# 

cp $HOME/lsrrinit fort.8 
# 

# RUN JOB 

# 

FILENV=.assign 

# Export the file environment: 
export FILENV 

# Assign statements for unformatted restart files coming from Iris: 
assign -a fort.8 -F f77 -N ieee -C ascii FORT8 

assign -a fort.9 -F f77 -N ieee -C ascii FORT9 

# SSD files 

# Export the SSD environment variable 
export SSD 

# See if there's enough space on SSD for all the files; if not, 

# use disk for all the files. 

if setf -n 34868b $SSD/dummyfil 
then 

rm $SSD/dummyfil 
FILDIR=$SSD 
else 

FILDIR=$TMPDIR 

fi 

# The number in the -n option is how many blocks of 512 WORDS will be 

# written to the file; e.g. if 2 arrays of dimension (512,3) will be 

# written to the file, then the assign statement would have "-n 6." 

# 

#******note - there are only files for 2 rows here!!********* 

# 

# DEFINE SSD FILES FOR Q ARRAYS - 1 FOR EACH ROW AND GRID TYPE 

# TOTAL SIZE OF EACH OF THESE WILL BE nAIRFOILS*5*nia*nja*nka 
assign -a $FILDIR/fort51 -n 1499 -s u FORT51 

assign -a $FILDIR/fort52 -n 1499 -s u FORT52 
assign -a $FILDIR/fort53 -n 1 -s u FORT53 
assign -a $FILDIR/fort54 -n 1 499 -s u FORT54 
assign -a $FILDIR/fort55 -n 1499 -s u FORT55 
assign -a $FILDIR/fort56 -n 1499 -s u FORT56 

# DEFINE SSD FILES FOR QOLD ARRAYS - 1 FOR EACH ROW AND GRID TYPE 

# TOTAL SIZE OF EACH OF THESE WILL BE nAIRFOILS*5*nia*nja*nka 
assign -a $FILDIR/fort61 -n 1499 -s u FORT61 

assign -a $FILDIR/fort62 -n 1 499 -s u FORT62 
assign -a $FILDIR/fort63 -n 1 -s u FORT63 


4.11 


assign -a $FILDIR/fort64 -n 1499 -s u FORT64 
assign -a $FILDIR/fort65 -n 1 499 -s u FORT65 
assign -a $FILDIR/fort66 -n 1499 -s u FORT66 

# DEFINE SSD FILES FOR AJA ARRAYS - 1 FOR EACH ROW AND GRID TYPE 

# TOTAL SIZE OF EACH OF THESE WILL BE nia*nja*KMX 
assign -a $FILDIR/fort71 -n 300 -s u FORT71 

assign -a $FILDIR/fort72 -n 300 -s u FORT72 
assign -a $FILDIR/fort73 -n 1 -s u FORT73 
assign -a $FILDIR/fort74 -n 300 -s u FORT74 
assign -a $FILDIR/fort75 -n 300 -s u FORT75 
assign -a $FILDIR/fort76 -n 80 -s u FORT76 

# DEFINE SSD FILES FOR PERMANENT XYZ ARRAYS - 1 FOR EA ROW AND GRID 
TYPE 

# TOTAL SIZE OF EACH OF THESE WILL BE 3*nia*nja*KMX 
assign -a $FILDIR/fort41 -n 899 -s u FORT41 

assign -a $FILDIR/fort42 -n 899 -s u FORT42 
assign -a $FILDIR/fort43 -n 1 -s u FORT43 
assign -a $FILDIR/fort44 -n 899 -s u FORT44 
assign -a $FILDIR/fort45 -n 899 -s u FORT45 
assign -a $FILDIR/fort46 -n 239 -s u FORT46 

# DEFINE SSD FILES FOR WORKING XYZ ARRAYS - 1 FOR EACH ROW AND GRID 
TYPE 

# TOTAL SIZE OF EACH OF THESE WILL BE 3‘nia*nja‘KMX 
assign -a $FILDIR/fort81 -n 899 *s u FORT81 

assign -a $FILDIR/fort82 -n 899 -s u FORT82 
assign -a $FILDIR/fort83 *n 1 -s u FORT83 
assign -a $FILDIR/fort84 -n 899 -s u FORT84 
assign -a $FILDIR/fort85 -n 899 -s u FORT85 
assign -a $FILDIR/fort86 -n 239 -s u FORT86 

# DEFINE SSD FILES FOR TMP ARRAYS - 1 FOR EACH ROW OUTER AND TIP 

# TOTAL SIZE OF EACH OF THESE WILL BE 5*nia*nja‘nka 
assign -a $FILDIR/fort91 -n 1499 -s u FORT91 

assign -a $FILDIR/fort92 -n 1 -su FORT92 
assign -a $FILDIR/fort93 -n 1 499 -s u FORT93 
assign -a $FILDIR/fort94 -n 1499 -s u FORT94 

# DEFINE SSD FILE FOR STR ARRAY - ONLY ROTOR INNER GRID 

# TOTAL SIZE WILL BE 5*nia*nja*nka 
assign -a $FILDIR/fort97 -n 1 499 -s u FORT97 

# DEFINE SSD FILE FOR INTERP ARRAYS - 1 FOR EACH ROW 

# TOTAL SIZE OF EACH WILL BE 3*(4*nka*nia*3 + nka*nia*3) 
assign -a $FILDIR/fort31 -n 435 -s u FORT31 

assign -a $FILDIR/fort32 -n 435 -s u FORT32 

# DEFINE SSD FILE FOR H-RATIO ARRAYS - 1 FOR EACH ROW 
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# TOTAL SIZE OF EACH WILL BE (nia*nja) 
assign -a $FILDIR/fort34 -n 1 -s u FORT34 
assign -a $FILDIR/fort35 -n 1 -s u FORT35 

# 

time ./a.slvr 
# 

# cat the unsteady Ps and Pt files 
cat fort.21 

cat fort.22 
cat fort.24 
cat fort.25 

# 

# COPY OUTPUT FILE TO CRAY 
cp fort9 $HOME/lsrr200 

Is -I $HOME 

# 

# SEND RESTART FILE TO FRONTEND 

# 

# dispose fort9 ????????? 


The dispose command must be edited in the file ’skeljob’ to send the output file back 
to the frontend. 

For two-dimensional cases in which H-ratio is to be specified, it is necessary to 
edit the here document defined as ’fort.4’ in the above jobstream. The format is the 
same as is already contained in the here document: for each airfoil, simply edit the 
number following ’nvals:’ to tell the code how many values of x and h to read (must 
be a two-digit number in columns 9 and 10, not greater than 25 for each airfoil), then 
list that many lines of format (FI 0.0) to define h as a function of x. The x-values 
should be in increasing order, beginning with a value that is not greater than the inlet 
boundary of the airfoil’s outer H-grid, and ending with a value that is not less than the 
exit boundary. The values of h can be defined by using results from either a 
streamline or a three-dimensional inviscid analysis of the airfoil row (Ref. 9). For 
two-dimensional cases with no stream tube contraction, no modifications to the here 
document are necessary, assuming the solution domain does not lie outside the range 
(-100.0<x> 100.0). 

In addition to the restart file, the output file from the UNICOS jobstream will be 
sent to the frontend. This file will contain all output from writes to Unit 6 in the 
solver, including convergence information and listings of the files containing midspan 
pressure distributions and exit total pressures. 

Experience with the solver has revealed several somewhat common problems 
that occur when a calculation is begun from an initial guess. The following is a list of 
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the most common, along with some possible causes and solutions. For any problem, 
the first thing that should be checked for is a grid error at the trouble spot, such as a 
discontinuity or lines crossing over each other. The solutions listed here assume the 
grid is not the problem. 


Symptom 

"TROUBLE AT GIVTRr 
message at beginning of mn 


Density and/or energy less than 
zero on or near the airfoil surface 
(particularly the trailing edge); 
may appear in the form of an error 
message from CONVRG, or a bad 
argument to mathlib occuring in MUKN 


Bad argument to mathlib occuring in 
MUTR routine 


Bad argument to mathlib occuring in 
CORREC 


Possible cause/cure . 

Outer H-grid lines defined by variables 
IBEG, IEND, JBEG, and JEND are not 
contained entirely within the inner 
O-grid; regenerate grid (In some cases, it 
is possible to visually determine what the 
values of IBEG, etc should be and simply 
hardcode them in the routine RESTIN.). 

Time step too large; reduce it and 
continue running. 

Incorrect Reynolds number has also been 
a cause; this can be hardcoded in 
routine RESTIN. 


Incorrect Reynolds number, or 
initial guess problems; less frequent 
eddy viscosity calculation in beginning 
has been known to smooth out "bad" 
points. 

Incorrect or discontinuous inlet 
boundary values; re-initialize flowfield. 
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RAI3DS Subroutines 


RAI3DS (main routine): 
RESTIN: 

JACBIN: 

TRIANG: 

GIVTRI: 

HDATA: 

EIGEN: 

CONTRL: 

GMOVE: 

RHSSI: 

MUKN: 

MUTRSI: 


MUTR2D: 

FLUXR: 

VFLUX: 

STCONT: 

LHSSI: 


reads input file in Unit 3; calls: 
reads initial flow file or restart file from previous run 
calculates Jacobians for the variable transformation, i.e. 
(x,y,z) -> (£,t|,<;) 

calculates factors for interpolation of flowfield at inner 
O-grid/outer H-grid interface; calls: 
for the given point in the given grid, find the triangle in the 
other grid in which the point lies and determine the weights 
of the three points for interpolation 
(called in 2D cases only) reads H-ratio data in Unit 4, 
linearly interpolates in x to define H at every grid point 
(called in steady-state cases only) calculates grid-dependent 
time step 

controlling routine for each solution iteration; calls: 

(called in multirow cases only) moves rotating grids relative 
to stationary grids according to size of time step; if rotors 
have passed stators, i.e., a cycle has been completed, 
applies periodic shift to rotors to begin new cycle 
controlling routine for computing right hand side of Eq. 3 
(chapter 4) for stator inner O-grids; loop 10 calculates 
^-direction contribution, loop 80 calculates T| -direction 
contribution, loop 160 calculates q-direction contribution; 

calls: 

calculates kinematic viscosity for stator inner O-grid 
(called in 3D cases only) calculates eddy viscosity for stator 
inner O-grid in q -direction and q-direction; loop 10 
calculates viscosity below midspan away from the wall; 
loop 210 calculates it above midpsan away from the wall; 
loop 410 calculates it below midspan close to wall; loop 
610, above midspan close to wall 
(called in 2D cases only) calculates eddy viscosity for stator 
inner O-grid in q -direction only 
computes the numerical flux vectors if, d , and (f in Eq. 3 
(chapter 4) using Roe’s scheme for stator inner O-grid 
computes viscous flux vectors S and T for stator inner grid 
(called in 2D cases only) calculates stream tube contraction 
terms (H-ratio) instead of the d vector 
controlling routine for computing the left hand side of Eq. 3 
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SMATRX: 

VMAT: 

BTRI: 


(chapter 4) and solving for Q ; loop 10 inverts the 
solution matrix in the ^-direction, loop 120 inverts in the 
T| -direction, loop 250 inverts in the ^-direction; implicit 
inlet/exit boundary conditions are applied within the 
^-direction loop, implicit surface boundary conditions are 
applied within the T|-di^ct^on loopj calls: 
computes the matrices A ,B , and C in Eq. (3) for stator 
inner grids; called once in each loop 
computes the matrices M and N in Eq. (3) for stator inner 
grids; called only in T|- and (^-direction loops 
inverts the given matrix; called once in each loop 


RHSRI: 

MUKN 

MUTRRI 

MUTR2D 

FLUXR 

VFLUX 

STCONT 


see above descriptions under RHSSI; all are the 
same here except they are applied to rotor inner 
O-grids instead of stator grids 


LHSRI 

SMATRX 

VMAT 

BTRI 


see above descriptions under LHSSI; all are the 
same here except they are appiled to rotor inner 
O-grids insead of stator grids 


PGRID: 


RHSSO: 


SRINT: 


MUKN: 

MUTRSO: 


FLUXR: 


defines solution above and below periodic boundaries of the 
outer H-grid, i.e. the boundaries in the circumferential 
direction, for use in RHSSO 

controlling routine for computing right hand side of Eq. 3 
(chapter 4) for stator outer H-grids; loop 10 calculates 
^-direction contribution, loop 180 calculates T| -direction 
contribution, loop 340 calculates ^-direction contribution: 
calls: 

interpolates flowfield onto outer H-grid at inner O-grid 
interface from the inner grid flowfield, using interpolation 
factors calculated in TRIANG 
calculates kinematic viscosity for stator outer H-grid 
(called in 3D cases only) calculates eddy viscosity for stator 
outer H-grid in ti - direction and ^-direction; loop 10 
calculates viscosity below midspan; loop 210 calculates it 
above midpsan 

computes the numerical flux vectors 
(chapter 4) using Roe’s scheme for stator outer H-grid 


T? , and Cj in Eq. 3 
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STCONT: 

VFLUX: 

LHSSO: 


SMATRX: 

VMAT: 


BTRI: 

PGRID: 

RHSRO 

SRINT 

MUKN 

MUTRRO 

FLUXR 

VFLUX 

STCONT 

LHSRO 

SMATRX 

VMAT 

BTRI 

RHSRL 

MUKN 

FLUXR 

VFLUX 

LHSRL 

SMATRX 

VMAT 

BTRI 


(called in 2D cases only) calculates stream tube contraction 
terms (H-ratio) instead of the (f vector 
computes viscous flux vectors and 'f for stator outer grid 
controlling routine for computing the left hand side of Eq. 3 
(chapter 4) and solving for Q ; loop 10 inverts the 
solution matrix in the ^-direction, loop 190 inverts in the 
q-direction, loop 330 inverts in the q-direction; implicit 
inlet/exit boundary conditions are applied within the 
^-direction loop, implicit surface boundary conditions are 
applied within the q-direction loop^ calls: 
computes the matrices A ,B , and C in Eq. (3) for stator 
outer grids; called once in each loop 
computes the matrices M and N in Eq. (3) for stator outer 
grids; called only in T|- and q-direction loops 
inverts the given matrix; called once in each loop 
see above description 

see above descriptions for RHSSO; all are the 
same here except they are applied to rotor outer 
► H-grids instead of stator grids 


see above descriptions for LHSSO; all are the 
same here except they are applied to rotor outer 
H-grids instead of stator grids 


see above descriptions for RHSSI and LHSSI; 
all are the same here except they are applied to 
tip clearance grids; note there is no call to an eddy 
viscosity routine - eddy viscosity is assumed to be 
zero in the clearance region 


CORREC: 


WALTMP: 

BCSURF: 


applies explicit boundary conditions at all boundaries, 
including all solid surfaces, periodic boundaries, inlet/exit 
boundaries, inner/tip grid interface for cases with tip 
clearance, and row interfaces for multirow cases; calls: 
calculates temperature on solid surfaces; called only if the 
option was chosen to specify wall temperature 
imposes solid surface boundary condition, either adiabatic 
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CONVRG: 

WRTPS: 

WRTEXT: 

RESTOT: 

PLOT3D: 


wall, specified wall temperature, or specified heat flux 
calculates and prints out change in solution from previous 
time step, specifically, change in energy 
prints out blade surface static pressure distribution at 
midspan 

prints out exit total pressure across pitch at midspan 

writes restart file with new solution 

writes grid and solution in PLOT3D format for viewing 


Miscellaneous routines: 


SHIFT: 


PUTXYZ,GETXYZ: 

PUTAJA.GETAJA: 

PUTQ,GETQ: 

PUTOLD.GETOLD: 

PUTH,GETH: 

PUTSTR,GETSTR: 

PUTTMP,GETTMP: 

PUTTRP.GETTRP: 


applies periodic shift to airfoil grids for multiblade cases; 
this saves on storage, as only one copy of the grid must be 
stored 

puts/gets grid arrays on SSD 

puts/gets Jacobian arrays on SSD 

puts/gets current Q arrays on SSD 

puts/gets differential (QD) arrays on SSD 

puts/gets H-ratio arrays on SSD 

puts/gets inner QD grid (for use in RHSRL) on SSD 

puts/gets Q or QD arrays to be stored temporarily on SSD 

puts/gets outer/inner interface interpolation factors on SSD 
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Symmetry Condition 


A symmetry condition may be used in the solver for axisymmetric cases in 
order to, in effect, increase grid resolution in the spanwise direction. There were two 
options for imposing symmetry implemented in the solver, only one of which was 
successfully used. The other was left in with the intention of later debugging, but the 
problem was never resolved and this option should not be used. This is symmetry 
option number 2 in the RAUOB shell script. The intent of this option is to allow the 
full grid capability of the program to be used for half the span, and modify the 
boundary condition applied at the uppermost K-line to assume a symmetry condition. 
That is, the geometry contained in the computational grid is only half of the actual 
airfoil. The difference between this and option number 1, which was run successfully, 
is that option 1 requires the full airfoil to be present and accounted for by several 
K-lines. In this case, most, but not all, of the K-lines are distributed below the 
midspan, the solution is calculated from hub to midspan, and a symmetry condition is 
applied at the midspan. 

To generate a grid to which the symmetry option is applicable, the user should 
have a three-dimensional axisymmetric configuration in the input file, and run the grid 
generator as usual up to this point: 

ENTER KMAX (RADIAL PNTS) BEFORE ADDITION, FORMAT=l2 
(DEFAULT=25, MAX=51(NO ADDITION)) 

The number of K-lines specified here should be based on the total number of lines 
desired for the half-span. The number that is chosen here will determine how many 
lines will be retained in the UPPER half of the span, that is, those that will not be used 
in the calculation. The remaining K-lines allowed in the solver will be distributed 
between the hub and midspan. For example, choosing 15 here will result in a final 
grid with 7 K-lines above midspan and 42 lines below midspan (due to the maximum 
K-dimension being set to 49 in the current version of the RAI3DS system). 

Continuing with this example, in which a three-dimensional airfoil with a total span of 
6.0 inches was used, 

15 

CURRENT DIMENSIONS: 

NOTE THAT POINTS WILL BE ADDED IF AXIAL SPACING AT 
INLET/EXIT BOUND IS CHANGED, AND AGAIN IN OVRLAP, 

ALSO RADIAL DISTRIBUTION CAN BE MODIFIED LATER. 

AIRFOIL: 1 

OUTER: 50, 31 INNER: 101, 21 RADIAL: 15 15 
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PAUSE- 


ENTER 1 TO INCREASE CLUSTERING AT T.E. FOR AIRFOIL 1 

ENTER 1 TO CHANGE K-LINE SPACING AND/OR DIST. IN TIP 
(DEFAULT: SETA1 = SETA2 = 0.03% SPAN = 1 .80000E-03) 

1 


This option must be excercised to re-distribute the K-planes. 

ENTER 1 TO DISTRIBUTE K-LINES FOR A SYMMETRY CONDITION 

1 


The above question will be asked only for a three-dimensional, planar stator. 


SPAN = 6.00000 KMAX = 15 

ENTER DIF BTWN HUB AND 2ND K-LINE, FI 0.0 

.002 

With this input, the program will use a spacing of 0.002 inches at the hub, and a 
spacing of 6.0/(NKA-KMAX) at the midspan, where NKA is currently set to 49, and 
KMAX for this case is 15. 

RADII OF K-PLANES ACCORDING TO ORIGINAL STRETCHING 
HAVE BEEN WRITTEN TO A FILE. 

HALT EXECUTION OR HIT ENTER TO CONTINUE 
PAUSE.. 


At this point, the plot of radii as a function of K-index looks like this: 




Note that there are two points above midspan that mirror the two points below 
midspan. This is necessary for the implementation of the line of symmetry boundary 
condition in the solver. The following is a view in the axial direction of the final grid 
generated by continuing the above example: 



Note that K-line 40 is at Z=3.0, which is the midspan. This is what should appear in 
the first line of FORTRAN unit 3, input to the flow solver. Unit 3 is defined by the 
here document ’fort.3’ in the UNICOS jobstream created by RAUOB. If the value of 
KMID (in this case, 40) is not properly specified, the symmetry condition will not be 
implemented properly. 
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Boundary Conditions 


Inlet Boundary Condition 

In the following expressions, u,v,w are velocities in the x-, y-, and z-directions, p is 
static pressure, p is density, c is speed of sound, s is entropy, and h is total enthalpy, 
a is arctan(vlu) and <p is arctan(wlu)', boldface type indicates a new updated value, 
subscript 0 indicates an initial value, and subscripts of 1 and 2 indicate a current value 
at inlet and one point downstream of inlet, respectively. 


The inlet boundary condition in ROTOR3 was the following: 

2 f'lyp} 

Reimann invariant 1 = R = u+ -h — °i 

1 0 Y-l J 


Reimann invariant 2 = R = u 

2 


2 ( 

2 " Y-l C^P J 


u x = (R + R 2 )/2 

V = V 
1 0 

w = w 


c = -^-(R -R) 
1 4 1 2 y 

s = s n 
1 0 


, C l/(7-l) 

P =( 4 “) 

1 y s l 


p c 

= r i i 


The following modifications were made to, in effect, hold inlet total pressure, in 
addition to reducing the stiffness of the boundary condition: 
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First define the variable h' as follows: 


h' = — + u = 2h - (v + w ) (1) 

7-1 

Assuming h = h Q , that is, inlet total enthalpy is held at its initial value, 

K =2 h - (v 2 + w 2 ) (2) 

1 0 1 1 


It can be shown that (2), together with the assumption 


u 

i 


2^ 


2c ? 

y- 



( 3 ) 


yields: 


_ 2 _ 

Y-l 


(14 



+ 


4R 

T 3- c 

y-l 1 


(h' -R 2 ) = 0 
1 2 


This is a quadratic equation for that can be solved with the quadratic formula. 
Then, 


v = u tan(<x ) 
1 10 

w = u tan{i p ) 
i i Mr 


s = s n 
1 0 


, C 1/(Y-1) 

P =( *“) 

1 


P i 


P c 

= r i i 
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Surface boundary conditions 


The surface boundary condition in ROTOR3 was a first-order adiabatic wall 
condition, that is, in addition to the no-slip condition for velocity, zero normal 
pressure and density gradients were imposed. Options were added to specify either 
wall temperature or heat flux on the surfaces, which modified the density condition. 

The specified wall temperature condition that was implemented was developed 
by Griffin (Ref. 7). It requires an additional input file containing adiabatic wall 
temperatures on the surfaces, and the wall temperature is prescribed as a percentage of 
the adiabatic wall temperature. The percentage is input in the RAIJOB shell script. 
The adiabatic wall temperature file can obtained by running RAI3DS with the 
adiabatic wall condition; the file will get written out in routine WALTMP, which will 
also read the file if the option to specify wall temperature is chosen. The UNICOS 
jobstream will need modification to assign these files to the file environment. 

The following is the derivation of the second-order density condition allowing 
specified heat flux. Here, Q is heat flux, k is the coefficient of thermal conductivity, R 
is the specific gas constant, T is temperature, p is static pressure, p is density, and y is 
distance from the wall. 

Assume 


Q = -k- 


dT_ 

dy 


a Jt 

dy 



Assuming a zero normal pressure gradient, this yields 

2 

kp dp dp pQ 

Q = — j-( — ) => — = , where Q = 

Rp dy dy p 

Now assume that p is some quadratic function of y, given by 
2 


p = ay + by + c 

Then, values of p at yl, y2, and y3 are given by the system 


QR 

k 
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2 

ay^ + by x + c = P 1 
2 

ay + by + c = p 
y 2 y 2 r 2 

2 

ay 3 + by 3 + c = p 3 

Assuming y is zero, that is, at the wall, the system reduces to 

^ 1^3 ^3 ^ 1^2 

D = 2 2 

V 3 - ^ 


Differentiating the quadratic equation for p, 


^- = lay + b 
dy 


From the previous expression for 


~r- = b at the wall. 

dy 

dp 


dy' 

(P 2 - p,)y 3 - (P 3 - p,)y 2 


V 3 


*2 >3 


p 2 Q' 


Algebraically solving for p^ 


p. = 


P/, 

y? - y? 


p Q’ (y ? 3 - y 2 y s ) 
P(y*-yp 


The values of p and p at y is used in the implementation of the condition. To define 
the value of Q' (note the solver will expect die given heat flux value to be properly 
nondimensionalized and to contain the factors of R and k), divide the dimensional Q 
by the dimensional k , and use a value of 1 .0 for R to be consistent with the solver’s 
nondimensionalization scheme. This will result in a value with units of degrees over 
length. Divide by inlet total temperature and convert remaining length unit to inches 
to obtain the value of Q' in inches . This is what should be input to the RAUOB shell 
script. 
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